{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook takes the heght anomalies computed using Arctic DEM and estimate  a linear gridded time series using a 3-D Guassian interpolation. This is modeled after the Science Data Generation tutuorial from ICESAT-2 Hackweek 2020."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "import numpy as np\n",
    "import cartopy.crs as ccrs\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from utils import transform_coord\n",
    "from utils import make_grid\n",
    "from utils import mad_std\n",
    "from utils import spatial_filter\n",
    "from utils import interp2d\n",
    "from utils import tiffread\n",
    "from utils import binning\n",
    "from scipy.ndimage.filters import generic_filter\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "usage: interpgaus3d.py [-h] [-b w e s n] [-d dx dy] [-t tmin tmax dt]\n",
      "                       [-r radius] [-a alpha_d alpha_t] [-p epsg_num]\n",
      "                       [-s n_sample] [-c dim thres max] [-v x y z t s]\n",
      "                       ifile [ifile ...] ofile [ofile ...]\n",
      "\n",
      "Spatio-temporal interpolation of irregular data\n",
      "\n",
      "positional arguments:\n",
      "  ifile               name of input file (h5-format)\n",
      "  ofile               name of ouput file (h5-format)\n",
      "\n",
      "optional arguments:\n",
      "  -h, --help          show this help message and exit\n",
      "  -b w e s n          bounding box for geograph. region (deg or m), optional\n",
      "  -d dx dy            spatial resolution for grid (deg or km)\n",
      "  -t tmin tmax dt     temporal resolution for grid (months)\n",
      "  -r radius           search radius (km)\n",
      "  -a alpha_d alpha_t  spatial and temporal corr. length (km and months)\n",
      "  -p epsg_num         EPSG proj number (AnIS=3031, GrIS=3413)\n",
      "  -s n_sample         sample every n:th point in dataset\n",
      "  -c dim thres max    dim. of filter in km, sigma thres and max-value\n",
      "  -v x y z t s        name of varibales in the HDF5-file\n"
     ]
    }
   ],
   "source": [
    "!python ./interpgaus3d.py -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, This runs a 3-D Gaussain interpolation with grid size and search raidus of 3 km with a time step of a half a week. These variables (as well as the model paramenters a and c), can and should be experimented with(we should probably increase the dt given teh amount of daat we have)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parameters:\n",
      "('ifile', ['/home/jovyan/shared/data-crossovers/Arctic_DEM_dh_filtered.h5'])\n",
      "('ofile', ['/home/jovyan/shared/data-crossovers/AD_filtered_data_cube.h5'])\n",
      "('bbox', [None])\n",
      "('dxy', [3.0, 3.0])\n",
      "('time', [2018.8, 2020.25, 0.125])\n",
      "('radius', [3.0])\n",
      "('alpha', [10.0, 0.125])\n",
      "('proj', ['32607'])\n",
      "('n_sample', [1])\n",
      "('filter', [100.0, 3.0, 10.0])\n",
      "('vnames', ['lon', 'lat', 'dh', 't_year', 'dummy'])\n",
      "-> reading data ...\n",
      "-> cleaning data ...\n",
      "-> creating kdtree ...\n",
      "-> interpolating data ...\n",
      "-> saving predictions to file...\n",
      "\u001b[0m"
     ]
    }
   ],
   "source": [
    "!python ./interpgaus3d.py /home/jovyan/shared/data-crossovers/Arctic_DEM_dh_filtered.h5 /home/jovyan/shared/data-crossovers/AD_filtered_data_cube.h5 -d 3 3 -t 2018.8 2020.25 0.125 \\\n",
    "-r 3 -a 10 0.125 -p 32607 -c 100 3 10 -v lon lat dh t_year dummy -s 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = h5py.File('/home/jovyan/shared/data-crossovers/Arctic_DEM_dh.h5', 'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X                        Dataset {17, 19}\n",
      "Y                        Dataset {17, 19}\n",
      "Z_nobs                   Dataset {141, 17, 19}\n",
      "Z_pred                   Dataset {141, 17, 19}\n",
      "Z_rmse                   Dataset {141, 17, 19}\n",
      "epsg                     Dataset {SCALAR}\n",
      "time                     Dataset {141}\n"
     ]
    }
   ],
   "source": [
    "!h5ls /home/jovyan/shared/data-crossovers/AD_filtered_data_cube.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#read in interpolated data\n",
    "with h5py.File('/home/jovyan/shared/data-crossovers/AD_filtered_data_cube.h5','r') as f_c:\n",
    "    Xi = f_c['X'][:]\n",
    "    Yi = f_c['Y'][:]\n",
    "    ti = f_c['time'][:]\n",
    "    Zi = f_c['Z_pred'][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "hi = np.nanmean(np.nanmean(Zi,1),1) #mean elevation chnage for each time step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Elevation residuals (m)')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtAAAAEGCAYAAABM2KIzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd5xU5fXH8c9hWemCUiy0RcSCYsVu/NkbGrux19iSWOIviUaN0SQm0STGGoldf9bYS0AxxhK7oAiKSFHpShMB6ezz++OZyb0zzO7O7M6de2fm+3695rV36n32ssycOfc85zHnHCIiIiIikp9WcQ9ARERERKScKIAWERERESmAAmgRERERkQIogBYRERERKYACaBERERGRArSOewCF6tatm6urq4t7GCIiIiJS4UaNGjXXOdc9+/ayC6Dr6uoYOXJk3MMQERERkQpnZlNy3a4SDhERERGRAiiAFhEREREpgAJoEREREZECKIAWERERESmAAmgRERERkQIogBYRERERKYACaBEpzMKF8PHH4FzcIxEREYmFAmgRyd+kSTBwIAwaBD/6kYJoERGpSpEF0GZ2t5nNNrOPm3jcDma22syOjmosIlIEX38NBxwAM2b460OHwh13xDsmERGRGESZgb4XOLCxB5hZDXAt8GKE4xCRllq0CA4+GD7/PPP288+H99+PZ0wiIiIxiSyAds69Dsxv4mHnA08As6Mah4i00IoVcNRR8MEH/nqrVtCvX+Z9c+fGNz4REZESi60G2sx6AkcAQ+Mag4g0ob4ezjgDXnopuO3vf/fXO3f216dNgxNOgNWr4xmjiIhIicU5ifAG4BLnXJOfumZ2tpmNNLORc+bMKcHQRASA666DBx8Mrl99Nfzwh9C/PzzwQHD7Sy/BVVeVfHgiIiJxMBfhLHozqwOed85tmeO+LwBLXe0GLAHOds493dhrDh482I0cObLIIxWRNXz0EeywA6xc6a+fcw7cdhuYBY+54gq45hq/3bo1fPUVdO1a+rGKiIhEwMxGOecGZ98eWwbaOdfPOVfnnKsDHgd+1FTwLCIlsnw5nHxyEDzvuCPccktm8Aw+I73ppn571SqYMKG04xQREYlB66he2MweBvYEupnZdODXQC2Ac051zyJJduWVMHas327XDu6/32eYs9XU+AD6s8/89XSLOxERkQoWWQDtnDu+gMeeFtU4RKRAb7wBf/pTcP2664Iscy69egXb06dHNy4REZGE0EqEIhJYtAhOOSVYYXDfff2Kg43p2TPYVgZaRESqgAJoEQn84hfwxRd+u3NnuOce3/e5MQqgRUSkyiiAFhFv1Cjf4znt1lszyzMaogBaRESqjAJoEfElGxdeGJRuDBniF0fJhwJoERGpMgqgRQQefRTefNNv19bC9dev2bKuIdkBdIS95UVERJJAAbRItfvuO/j5z4PrF14Im2yS//PXXhs6dvTby5bBN98Ud3wiIiIJowBapNpdd13Qfq5HD7+6YKFUxiEiIlVEAbRINZsyxQfQab//ve++USgF0CIiUkUUQItUs1/8wpddAGy/PZx+evNeRwG0iIhUEQXQItXq9dfhH/8Irt94Y9M9nxsSDqC1GqGIiFQ4BdAi1Wj1aj9ZMO3442G33Zr/espAi4hIFVEALVKN7r4bRo/22+3awbXXtuz1wguuKIAWEak+M2bAb38Lb70V90hKQgG0SLVZsAAuuyy4fuml0Lt3y15TGWgRkeq1ZAnsvTdceSXsvrufnF7hawIogBapNr/5Dcyd67f79IGf/azlr6kAWkSkel1+OUyY4Ledg0sugZNPhqVL4x1XhBRAi1SKuXPh448b/9Y/fjzcfHNw/U9/gvbtW77v9daDmppgHMuXt/w1RUQk+f7zHz8JPduDD8Iee1RsUkUBtEglmDfPrx44aBD85S8NP+7ii2HVKr+9xx5wzDHF2X9NDay/fnB95szivK6IiCTXd9/59qfpxM0BB8DZZwf3jxwJW2wBRxzhP5vefRdWrIhnrEWmAFqkEvz738ES2lddFZRohA0bBsOH+20znzEwK94YVMYhIlJdLrsMJk/22507w513wtChcMstwVnJb7+Fp5/25YI77wxdusCQIf4xn38e39hbSAG0SCUIB6zffQfXX595/4oV8NOfBtfPOgu22aa4Y1AALSJSPV57DW66Kbh+ww2+I5MZ/PjHMGKEn2eTbelSn9A5/3zo3x823RQuughefDFY2KsMKIAWqQTZAevNN/uyjrRbbgkmeHTuDL/7XfHHoABaRKQ6rFgBZ5wRXB8yBE49NfMxe+8NX3wBn3wCf/+7n1RYV7fma02Y4M+IHnggrLuuf62bb4YxY6C+PtJfoyVaxz0AESmC7NX/Fi/2WehrroGvv4arrw7u+/WvoXv34o9BqxGKiFSHt98Oyi+6dPEBcq6SwFatYOBAf0nXRn/+uS8nHD7clx+GO3Wks9PDhvnr66wD3/uen7Ozxx6w3XZBaUjMlIEWqQS5Mr433eSz0FdcAQsX+ts23dSfWouCMtAiItVh/Phg+/vfz3z/b8pGG/nPoeefh/nzfenGRRf5z6ds33wDzz4b1E8vXtzysReJAmiRShAOWLt29T8XL4Yf/hDuuiu4769/hbXWimYMCqBFRKrDZ58F25ts0vzXadsW9t/ffzaNH+8nJN56q+8Qtd56mY/dZhtfgpgQKuEQKXfOZQas117rA2fwM5/TDj4YDjoounEogG7aypVQWxv3KNZUXw/vvOPrFefODS6NLYLQUAeXvn39Igrt2kUzVhGJX3pODeTOHDfXRhvBj37kL875/bz+ur8MHFi8/RSBAmiRcjdvXrBwydprw2mn+frnceOCx7Ru7b/hRykcQM+c6d/8itkmr9z96lfw+9/7iTd33BH3aDL9/Odrdm5pia++8q2sRKQyFSsD3RgzH5xvuqnvHJUwKuEQKXfhbG/Pnn6Cxa9+lfmYCy+M7k0urWPH4PTaihW5e1FXq8mT/YTO+nrfJzVJx2bRIrjttuK+5h13wNixxX1NEUmGFSv82aq0AQPiG0uMlIEWKXfhALpXL//zmGP8qk8jR/rbsgPqqPTs6Zvmp8cVRbePcnTTTZlLrE+ZAt26xTeesCefDEo1evb0K4Z16+Zr6Tt0KOwswt13+2V96+t93/GXXtJZCJFK8/nnsHq13+7Tp2rLtRRAi5S7cMu4dBlFTY1vYv/CC7DnnqWbeNGzZ1A6MmNG8RdrKUcLFvjAMmzqVNh++3jGk+3++4Ptiy7ys92ba4cdYOut/Yfryy/7WfaHHtryMYpIckRV/1xmVMIhUu6ySzjS1lkHjj8eNtigdGPRRMI13Xnnmq2Xpk6NZyzZpk6FV17x261awQkntOz1ttgi6PUKPhhfsaJlrykiyVKK+ucyoABapNw1FEDHQQF0ppUrM5e6TUtKAP3gg0Fpyb77woYbtvw1r746OOMxYQL87W8tf00RSQ5loIEIA2gzu9vMZpvZxw3cf5iZjTGz0WY20sx2j2osIhUtXMKRroGOi1YjzPTEEzBt2pq3JyGAdg7+7/+C66ecUpzX7d49s+b+6qszl5UXkfIWzkArgI7EvcCBjdz/MrC1c24b4AzgzgjHIlK5lIFOJuf8RM60ffYJtpMQQI8aBZ9+6rc7dIDDDy/ea//kJ9C/v99esCBzKXkRKW8q4QAiDKCdc68D8xu5f7Fz/52W3gFwDT1WRBqhADqZ3nzTd0EBaNMG/vCH4L4kBNDhyYNHH+2D6GJp0wb+/Ofg+sMPF++1RSQ+CxbA7Nl+u00b34WjSsVaA21mR5jZeOCf+Cx0Q487O1XmMXLOnDmlG6BI0i1ZAt9847dra+NvG6cAOhBemOSkk2Dbbf1EPfALjaQXv4nDypWZQW2xyjfCDj00aGE3dy6sWlX8fYhIaYXrnwcMCN7TqlCsv7lz7inn3GbA4cBvG3nc7c65wc65wd3jDhBEkiQcpG64YfxvZj16+FUPwQf2jS0FXclGjMhcRv2nP/XHJTxJL84a8RdeCBZz6d3btzostpoa3wkmbX4DJyRfe80H8wqwRZJPEwj/K69PWzNrZWbbmtkQM9vbzNYr5iBS5R79zSwhKwuIlIkklW+AD+DDbfOqMQud7n2crlA78EDf3g0yT3fGWcYRLt848cTovnh17Rps55pIOHq0D95POCGz5ENEkkn1z//V6LummfU3s9uBScAfgeOBHwEvmdk7Zna6mTXrndfMNjbz5/fMbDtgLUBTtUUKkbQAGqq7jOOpp+DII4Pex717Z7ZxCwfQubpzlMKiRfDss8H1k0+Obl9NBdBvvRVshzuCiEgyKQP9X02tRPg74DbgnNCEPwDMrAdwAnAycF/2E83sYWBPoJuZTQd+DdQCOOeGAkcBp5jZSmAp8IPsfYhIE3KtQhi3cCu9agqgH33UZ3PTS9z26wf//jfU1QWPSUIGetKkIMDfZBMYODC6fTUVQKfLSMCvYPnll5nHS0SSRRno/2o0gHbOHd/IfbOBG5rz3NT91wLXNjVAEWlEOECNuwd0WjVmoO+4A849F+rr/fUBA/xS1r17Zz4uCQF0OJAtxsIpjQkH0OFguaHbhg+H886Ldkwi0jz19TBxYnBdGeimmVkNMASoCz/HOXd9Q88RkRJIYglHOGj69tv4xlEK9fVw+eXwxz8Gt22+uQ+ecy2hnoQAOhy0dot42kkhGWiAYcMUQIsk1YwZvvMT+P/b664b73hillcADTwHLAPGAvXRDUdECpLEALp9+2D7u+/iG0fUli2D007zpRtp227rO1z06JH7OUkIoMOBbDjAjUJTAXR2W9KXX/bHtW3baMclIoVT/XOGfAPoXs65rSIdiYgULknLeKeFA+h0tqLSzJsHhx3mF0tJGzIEHnkEOnZs+HnZAbRzQa/kUklyBnrpUt/W7oADoh2XiBRO9c8Z8u2gMdzM9o90JCJSmNWr/YIcaVHXs+YrvKJdJQbQn3wCO+2UGTyfd57v+9xY8AzQpUvwmO++CxbBKaUkZaBz1UUPGxbdeESk+ZSBzpBvAP0O8JSZLTWzhWa2yMwWRjkwEWnC118HHR+6dfPLqiZBJWegn34adt4ZJk/21818/+Jbbw0WkGmMWfxlHEkJoJ1TAC0ShQkT/GJO9UWuuFUGOkO+AfRfgF2A9s65tZ1znZxza0c4LhFpShLrn6EyA+j6erj6ajjiCFi82N/WoQM8/jj87/8WVoYRdwCdlBKOJUt8vTPAWmsFZy4mTcqc6S8i+fviCxg0yJdBXX11cV9bGegM+QbQE4GP1adZJEGSWP8MlRdAz5njF0e56qrgtn794O23/e2FijuATkoGOjyBsEcP2Hff4Lqy0CLN88wzQZ/3P/0JZs8uzusuX+77tINPGPTvX5zXLWP5BtCzgFfN7JdmdnH6EuXARKQJykBH77HH/DLczzwT3LbPPvD++z7L0xxxB9BxZqDDOZjscRx8cHBdAbRI84TnZixdCn/5S3Fed/LkoCSkrk6dcsg/gP4CeBm/3Han0EVE4qIAOjpffw1HHw3HHpuZKb3oIt+mriWZ2/DiKpWegW7Xzl8AVq3yy4inZQfQBx0UXH/11cpugSgSBefgrbcyb7v11txzDQql+uc15NXGzjlX5EIaEWmxJC7jDeXdB9o5eOghuOACmD8/uL1XL7/a4IEHtnwf4Qz0tGktf71CLF8e1HDX1EDnztHvs2vX4G913jxYOzV9JjuA7t0bttwSPv7Yn4L+97/h0EOjH59IpZgyBWbOzLztu+/gr3+Fa65p2WuPHBlsq/4ZaCIDbWa3m1nO85Rm1sHMzjCzE6MZmog0KonLeEP5ZqBnzvS9nU86KTN4PussH9QVI3iGeEs4srPPpehB3VAddDiA7t7d/1QZh0jzhbPP4S/HN9/cspaZTz6Zudrq1ls3/7UqSFMlHH8DfmVmn5rZY2b2NzO728z+A7yFL+N4PPJRSjLdcYfP1M2aFfdIqlNSSzjKrQ+0c3Dvvb7W+bnngtv79oWXXoLbby9uprZnzyBwnTkTVq4s3ms3pZTlG7n2E95/uDQmXYsdDqCHD492XCKVJlz/fMEFsPnmfnvRIrjxxua95osvwnHHBfXPW23lr0vjJRzOudHAsWbWERgMbAAsBT51zn3W2HOlwr3xBpx9tt9+7TV4552g1lGi51xyA+jw38GSJfGstleIs8+GO+/MvO1HP/IZl04RTPVo0wbWX99/8ayv90F0377F308upZxAmJZPBjo9ll13hVat/HGZMsWXcqy1VmnGKVLuwhno733Pl1qcdJK/fsMN8NOfFpYMeOMN37oz/SV/wADfXzp8lrGK5VsDvRh4NdqhSFl5+ulge8wYuPBCn6mT0vj226C+uF07v8JdUrRu7YOeFSt88Lx8eXJnbC9ZAnfdFVzfaCN/fc89o91vnz7BmZupU0sXQCcpA50rgK6thXXXDe6bP99/2RCRxi1a5D+LwX8J3WknH+hefbXvq/7tt3Dyyb57UKtWPqnRqlXmtnM+WF6xwl/uvNN38gD/nvWvf8F668X3OyZMXgG0yBpefDHz+h13wP/8D5yokviSyK5/TlqGt337oBfpkiXJDaBnzgxaq62/vv8ACpegRKVPH3j3Xb9dyjrocACbxAx0+vHp++bNUwAtko933w3KLAYNCibrXn45nHaa337uucwytXytt54PnsPzNyTvNnYigRkz/KSqbOecA+PHl3481Sip5Rtp5TKR8Kuvgu2+fUsTPEN8EwnDQWuSMtDpSYSQGUwXo/2WSDUI1z/vtluwfeKJQS10c6yzjp8LMmBA81+jQhWcgTazVkBH59zCCMYj5WDEiGB7l138h+KECb6k4Jhj/Ddh1UhFK6kt7NLKJYAOT4DdYIPS7TeuADpJJRy5JhE29ngRaVi4/nnXXYPt1q199vixx/xndH29P+tWX7/mtpkvv6ut9T/btYMhQ/zCKbKGvAJoM3sIOBdYDYwCOpvZ9c65P0U5OEmoF14Ito84Ag44wNdbLVvmM9N9+/rTR2ut5SdM1dT4x5oFl6au19X5RSwOPFCTiHIppwx0kntBhwPoUpYKJCEDHWcJR319w8G8AmiRwqxeDW+/HVwPZ6ABNtzQz1OSoso3Az3QObcw1fN5GHAJPpBWAF1tVq/2p3PSDjjAt7W5+WbfLxf8h3RLT72++SY8+KCfUHTssX7yQ/hbdbX79NNgO0k9oNOUgW5ctWegFywI6jXTX7bTVMIhUphPPglW+dxgg9JNSq5y+dZA15pZLXA48IxzbiXgohuWJNbIkUFD9g028JMVAM48E3784+Lvb/58GDrUf6P+wx+K//rlaNQoeOSR4Po228Q3loaUSy/ocA10tQXQcWagG8uEKwMtUpjs+uekTSqvUPlmoP8OfAl8BLxuZn0B1UBXo3D3jf33zyy/uOUW+N3v/DfhFSt8+7Lly4M6q/QFMq9n37Zihc9yP/hg5lLH//d/8Mtflub3jNLKlfD55z5QKHQ1uPp6/0UlfcwOPBB23z2acbaEMtCN69rV1xcuXQoLF/oWU6VYVjspkwgbmkDY0ONFpGHhAFpnaksm3z7QNwE3hW6aYmZ7RTMkSbRwAH3AAWve36VLcXoS7703XHMNvPoq7LOPv23yZF9Ckq6pLjdffQV//7vPqKczn2uv7XsP9+8PAwfC4MH+suGGuV/jrruC9mdrrQU33ZTMbEM5BtClrIE281noz1LrUU2dGpzNiVIcJRxdugSLo6S/XDc0gTD7uko4RJoWnkCYXf8skWk0gDazi5t4/vVFHIsk3YIFQfBmBvvtF+3+WrXygfR668HXX/sP3unTy6++a/Ro+NOf/Czo7GWbFy70948eDU88Edy+/vo+k/CLX/gJmuCDn0svDR5zySXJbS1UjgF0KTPQAL17BwH0tGnRB9CrVvn/w+D//66zTrT7S2vVyu8rHbzPn68SDpFimTULvvjCb7dtm8ySvgrVVA10pyYuUk1eftlngAG23750NZQbbxxsT5xYmn0Wy9ixsMMO8NBDmcHzuus23nP4q6/gySdh5519H8+pU335yvz5/v66umSXs5RDAL1yZRDImZV+ha3wF8HXXot+f+m/HfABbSnP5GQHxQqgRYojnH3ecUd1rSqhRjPQzrmrSzUQKQPh9nW5yjeiMmBAUOM1cSLsu2/p9t1Sw4f7zF/arrvCBRfAkUf6/pxz5vjSlIkTfRZ65Ej44IPM1m8PPeSD6eXLg9tuusnX0CZVOQTQs2cHteTdu/t/j1IaMiRYRvy22/zZhSizwnFMIExrLIDOroFWCUd+3n/fn7UKv7/ko00bf2Zvzz3LtxxOAg31f5bI5dsHui1wJrAF8N81eZ1zZ0Q0Lkka5zLrnw88sHT7DmegJ00q3X6L4ZNPgu0//tGXXYT16OEvu+wCp5zib1u92meuf/tbHziD77Gdduih/pJk5dAHOq7657TDDvMrhH36qa8NvvlmuPLK6PYXxwTCXPtrKgO97rrB9jfflPe8h6h89BHssUfm+0Ihfv973z/+hBN8i9BS1N9LNCZMCLa33z6+cVShfNvY/R+wPnAA8BrQC1gU1aAkgcaPDzpirL12UJdbCuE633Ir4QgH0DvvnN9zamp8HdsTT/hT++E3xbZt4cYbizvGKJRDG7s465/B1wZfdllw/cYbYfHi6PaXpAx0Y5MIW7cOJiI7F7TNFG/hQr/ia3OD57QZM/zcjK228vNZwme4pHx8/XWw3dDkc4lEvucsN3bOHWNmhznn7kutTPhik8+SyhHOPu+zj1/qs1TKtQa6vh7GjQuub7FF4a+xxx7w3nvw8MO+I8npp0O/fkUbYmTKoYQjrh7QYccdB7/+tW9rmO55/rOfRbOvODPQ4SC5qQw0+PGlJzzOm1f6gD+pnIMf/jB4H+zQAa64orDyoy+/hH/8I/NLzL/+Bc884xetkvIye3awXep5HFUu3/916dlPC8xsS+AroK6xJ5jZ3cAhwGzn3JY57j8Rv6IhwGLgPOfcR3mOR0otXGe1//6l3Xc4A/355/Gd0q2v972uFyyAn/+86RrkL77wfX7Bv7E1Nwho1cpPJDzxxOY9Pw7lEEDHnYEGH/hccgmcc46//pe/wE9+4s80FFscLexy7S/fAHry5ODx4t16q+/mk3bHHXD88YW/zl//CiNGwHXXweuv+9tGjFAAXW6cy8xA9+gR31iqUL4lHLeb2TrAr4BngXHAdU08516gsULZL4D/cc5tBfwWuD3PsUgcwrXHW29d2n136hR8s16xInNxlVJ6/nm48EKfMbzhhqYfHy7faE72uZyVWwAdRw102qmn+npU8Fnxu++OZj9JKuFobBIhaCJhLu+9BxeHOsued17zgmfwZxCHDMlc3XXEiGBSrZSHRYuCUp527aBjx3jHU2XyCqCdc3c6575xzr3mnNvIOdfDOTe0iee8Dsxv5P63nHPp4rZ38HXVkkTOZQbQ/fuXfgxJmEj4zjvB9ksvNf34jz8OthVAJ08SMtDguyKEyzauvXbNfuHFkJRJhF995VdeBH92JdfCS2pll+nbb312OP13sf32PovcUjvu6Oe0gE9MpPuSS3kIl2/06JHMRbUqWL5dOHJODXfO/aZI4zgTGN7I/s8Gzgbo06dPkXYpeZszx3/TBZ8NzpUxiloSWtmlm9UDjBrlSzpaNfIdNJyB3nKNKqbKVg4BdBJqoNPOOsuvvDl3ru/5vdFGvp+rWX6XNm38wix9+/pLv35+Ylh4MmdSMtDhrgFdu+b+P6QAOtNjj8GUKX67c2dfw9ymTctft3VrP6flqaf89REjYLPNWv66Uhrh8g3VP5dcvjXQ4T5UbfG1zZ8WYwCpJcHPBHZv6DHOudtJlXgMHjxY55hKLZzx3XjjeL7lJmEiYTiAXrjQj2PTTRt+vEo4vKQG0EnJQIMPdH/6U7j8cn99+vTCX2PUqMzrAwb4dojpQCspGejPPw+2GwrkVcKRacaMYPu88/wXrGLZf//MAPqCC4r32hItTSCMVb4lHH8JXa4B9gR6tnTnZrYVcCdwmHNOaYakSk/mgcxAtpTCEwnjKuEIB9DgFzJoyKpVvr9vmgLoZHEuMwMdZw102k9+UtwzFRMnZv6NJmUSYX19sN3Q2SxloDOF2xsWe7Gd8KTwV15RO7tyogmEsWru0lvtgRZ9BTazPsCTwMnOuQlNPV5ilJ2BjkPcvaAXL878tg9+1cCTTsr9+MmT/YRH8BPEctV5VrJw6UASF1KZPz/49+nUqfFl1Utl7bX9AhlTp/pOM87lf/nuO/+8KVPgkUfgww/9a37wAeyeOrmXlBKOsIbGoQA606LQsgudOhX3tTfayM9rmTzZf9l96y3Ya6/i7kOioQx0rPKtgR4LpEsnaoDuQKP1z2b2MD5T3c3MpgO/BmoBUhMQrwS6An8zXxKwyjk3uPBfQSIX9wRCyAzc42hl9+WXa97WWAa6mss3IPkZ6CTVP4e1agV1dc177i67+J9t2wYBdLqso77ef2lIC6/2Vwpt2vgvKdlfplTCkZ9wBjqKTgv77++XkwdfxqEAujwoAx2rfDPQh4S2VwFfO+dWNfYE51yj/XWccz8Efpjn/iVOSchAp1vZff110MquuYFGc2SXb4APUlatyr2IQTV34IDkB9BJqn8utu22C7Y/+MD/XLAgKJ3o3Lm0CyGlde2afwCtDHSmKDPQsGYAHW5vJ8mlSYSxarQG2szWNbN18ct2py9LgbVTt0s1SEIAnb3vUpdx5Aqgly7NXGkwrJo7cEB5BdBJqH8upq23Dib6jhvnj3+cEwgb268C6PxEnYHea6/gjN4HH2SuUijJpRKOWDU1iXAUMDL1cw4wAZiY2h7VyPOkUnzzTXDqt127eLN1cU4kzBVAQ8NlHNVewlFbG3wgr1oVTV/jlqjkDHSnTrDJJn67vt534ohzAmFj+813EmG1L/ARdQa6c+egBAj80t6SfCrhiFWjAbRzrp9zbiPgReBQ51w351xXfEnHk6UYoMQs3IGjf//G+x5HLc6JhOHWW9tsE2yPHLnmY1esyFyQYODA6MaVVGbJzkIntQa6WLLLOOKcQJhWSAa6bWWQRr4AACAASURBVNtgYueqVb5tZDWLOgMNmd04RoyIZh9SXMpAxyrfaGgH59yw9BXn3HDgf6IZkiRKEiYQpsW5GmE4A33MMcF2rgz0xIn+Qx/8ohZRZIzKQZID6ErOQMOaAXS5lXBk31ftEwnDGehSBdDVnvVPuuXL/dwG8Gf7Sj0xWPIOoOea2RVmVmdmfc3sckCFadUgKfXPEF8G2rnMAProo4PtMWPW7Jta7eUbaeUSQFdaDTT4pZ7TyjEDnf34aq+DDmego/pCPnhw0G5z5szM9zFJnnCdevfu8Z4drlL5HvHj8a3rngKeBnqkbpNKl6QAOlcru1KYNy/4AOvY0Qfy6bGsXOmD6LBq78CRluRe0JWegd5222B77FgfEKWVSwZaAbTnXGlKOGpqYN99g+sq40g21T/HLt+VCOc75y50zm2bulzonJvf9DOl7CUpgE63soOglV0phLPP/fr5+t7BoZbl2WUcykB7Sc5AV3oNdJcuwXLPK1fC668H9yUlAx2uc85FJRze8uVBSdhaa/lLVPbbL9gO/81I8qiFXeyaamN3Q+rnc2b2bPalNEOUWCUpgIZ4yjjCEwj79fM/d9ghuC17ImG1t7BLS2oAvWRJMCmttrZyawfDddCjQk2TkpKB7tYtaLfX1OOrOQNdiuxzWrgTxyg12ko0TSCMXVMLqfxf6uefox6IJNDixcG33Npa6N073vGAD+LfeMNvT5qUmTGJSnYGGjID6HAGetmyILA3g803j358SZXUADq7/rmxIK6cbbcdPP643w5PCEtSAJ3v46s5gI66hV3Y5pv7dqVLl8L06f79X8FZMqmEI3ZNtbEblfr5WvoCjAG+SW1LJQu3sOvXr7RLZzckjgx0rgB6222DSRvjxgU1vp99Fqz41q9fZhBZbcohgK7E8o20cAY6LCklHE2NQyUcXikz0K1b+4V40pSFTi5loGOXVw20mb1qZunVBz8C7jGz66MdmsQuaeUbEM9qhOEAOl1X2rFjkF2ur/fLeoPKN8KSGkBXev1zWngiYVhSMtANLaKS6/HKQHulaIkZnt+hADq5lIGOXb5dODo75xYCRwL3OOe2B/Zt4jlS7pIYQMexGmGuDDSsWcYxfTo8/3xwWzVPIITkBtDVkoHu0QN69Vrz9rgC6M6dM89iqYQjP6XMQENmC0QF0MmlDHTs8g2gW5vZBsCxwPNNPVgqRBID6FK3slu9GqZMCa7X1QXb4UzNZZf5GvGHHw5uUwAdbCuAjkd2GUeHDr77RRzMMidsqoQjP6VYRCUsHEDnWmlVkkEZ6NjlG0D/Br+c92Tn3PtmthFQ4rWUpeTCNdBJCaCzW9lNmBDt/mbM8G3AwJ9yDn+AhTPQy5ZlPq9DB9hnn2jHlnRJ7QNd6YuohGUH0HFln3PtXxno/JRiEZWw9ERC8O9/4UBNkkNt7GKXbx/ox5xzWznnzktd/9w5d1S0Q5PYJWkZ77Bwq6V//jPafTVUvgF+sk34jatNG78QwR//6GuhKz04a0pSM9DVUgMNmdlEiG8CYVohAXT4/nnzqndp6VKXcLRuDdtsE1xXGUfy1NdnrkSoDHQs8p1EuImZvWxmH6eub2VmV0Q7NInV0qXBQiWtWmWWLsTt+98Ptp+NuB15rgmEaW3awKuvwvXXw0svwTff+J+XXAJ9+0Y7rnKQ1AC6mks44s5A77ij/9mq1ZrBfbb27f3/MfBneJL0N1RKpZ5ECCrjSLp584JuT126RLu4jjSoqT7QaXcAPwf+DuCcG2NmDwG/i2pgErNw4Ni3b7L+gw4Z4uspnYM33/T1kVFl1hrLQANstpm/yJoUQMdvgw38WZL06d64M9BXXuknNm65ZdNntcx8wJ9ehnzevMZXLqxUpc5AgzpxJJ0mECZCvjXQ7Z1z72XdtqrYg5EESeIEwrQePWDXXf12fX20ZRy5ViGU/CQxgF61Kjj1aVb5pz7NMrPQcWegu3SBiy+G/ffP7/HZZRzVKO4MtALo5NEEwkTIN4Cea2b9AQdgZkcDsxp/ipS1JAfQULoyjqYy0NKwJAbQX38d1NJ26+ZX2Kx06bIJgD594htHc4QD/mrtxBFHBnqzzTInEobnDUj8lIFOhHwD6B/jyzc2M7MZwEXAuZGNSuIX7sCRpAmEaeEA+sUX1+yCUSwKoJsviQF0NU0gTDv/fN8RZp994LTT4h5NYdSJI54MdOvWmQvxKAudLOrAkQhN1kCbWQ1wnnNuXzPrALRyzi1q6nlS5pKegd50U7+oysSJvkXav/8NBx9c3H0sWxbUX7ZqVX7Zu7iF61WTEkBXU/1zWteu8K9/xT2K5lEv6Hgy0ODLON56y2+PGuXnnpSDWbNg6tTc9zXWyaVPH9hww2jGVGzhDLRKOGLTZADtnFttZtunthPUzFUiU18Pn30WXE9iAG3ms9B/+Yu//uyzxQ+gwwuo9O5dHaf7iymcgU5KH2id+iwvykDHk4GGxjtxzJzpz+asXLnmZcWK/G/v0MG/b2+1lX9Pb4kvv/STVB94oHktD83g6KPh8st9i9IkUwY6EfLtwvGhmT0LPAb895PQOfdkJKOS+NTXw49/HASPNTVrtm9LisMOCwLo556Dv/3NZ4qLRRMIWyaJJRwLFwbbXbrENw7JjwLo+DLQDXXiuPFG+OlPi9eX+7LL/NnEY46BY48tPJieMwd+9zu47bZg0avmcA4ee8xfDj0Urrgic/5AkmgSYSLkG0CvC8wD9g7d5gAF0JWkvh7OOw9uvz247bzzgskkSbPLLv4Ddt48nxH54IPMN/2WUv1zyyQxgP7222B77bXjG4fkRyUcpV/KO22zzfz/4SVL/PvrrFm+z/1FFxV/XxMnwu9/7y/HHAOPPppfEP2Pf8CZZ2Z+yQAYNKjhz61cr7t8OYweHVx/7jl/+cMf4NJL8/89SkVn0hIhrwDaOXd61AORmNXXwznnwJ13BredeCL89a/xjakprVv7urz77/fXn3lGAXSSJDGADmegFUAnnzLQpV/KO62mxq9ImK6DvuYaGDo0uL9nTz+PoLbWrxNQWxtcsq83dNv48b78Lvw7PvaYz3CHV5zNpb4eLrww87m77ALXXgvf+17hv++HH/oA/oknguz6bbclM4BWBjoR8s1ASyVatMgvOT12LDz/fGY7uJNPhnvu8W+iSfb97wcB9LPPwm9/2/LXdA7eeCOzv3RSy1iSLOkBdOfO8Y1D8tOSAPrLL/2XpHXXLeqQSi6uEg7wCYl0AH3rrcHtgwbB668Xpwxq6VLfSem66+Dtt/1tQ4c2HUB/+GHQVadzZ7jvPv950Nxa6m239cH7J5/4hX4Apk/3ZSFJmv/inDLQCaEAutKsWOFX5/viC1i92i8csXp10FFixgz/pjBtWrBUd7ZTT4W77kp+8AxwwAE+s7FiBYwZ42uiN97Yr3a2/voN/w4Nvcm+9x788Y/Bh0bappsWd9zVoG3bYHv5cv93GPfflEo4yktzSjg++MDXrw4f7gPO0aOT2YozH/X1mRNwS70SY67l1uvq4IUXijeHoF07OPxw/569ww7+tkcf9Wc/G/vyM3x4sH3IIX5OTDFssYXvxjFzpj/+06cn6wzk4sX+Swf4Y1fqL1XyXwqgK8GMGTBsmL/8619r1oMV4owzfA103IFOvjp29P1t02+mP/tZcV/fDH74w9wfJNI4s6CGEvybftxv9irhKC+FZKDHj/ddGB57LLht8WI/uTg92bjcZAfPxZwknY/s973u3WHEiGjavQ0e7Pc3apT/wn3ffb6UoyHhAPqgg4o7lrq6oIXpl18mK4DObmHX0u4l0mx5/W80szZmdoKZXWZmV6YvUQ9O8vDgg9C3L5x9Njz9dGHBc+vW/tv2ccf5WcyvvOJroMsleE47PYIS/dpaPznl00/9Fwq9STVP0npBq4SjvHTuHLwfLV7sA6tss2bBWWf597Jw8Jz28MP+7Ec5iquFXdpmm/nPF/BffocP9x0zonJuaH22oUMb7vQxfz68847fNvNnIosp/TuDD6CTRC3sEiPfDPQzwLfAKCDHO9iazOxu4BBgtnNuyxz3bwbcA2wHXO6c+3OeY5GwG25Y88OhXz/YdVd/eqemxgfKtbW+pKFXLz/5o1cv3zh+rbXiGXcxHX20P6X43nu+LGX6dH+ZMyf3G3Bj7ZfatYOjjoKLL/bHSVomaXXQKuEoL2Y+C53Ous2bF2Q/v/sO/vxnXzub/bd15JF+HsPs2T7AfuUV2Hff0o69GOKsfwb/+fHcc/D44z7Rsvnm0e7vuOPgf//Xf9GdMMH/u+2995qPe+klX14BvtVcuNSnGOrqgu0kB9CaQBirfAPoXs65Awt87XuBW4D7G7h/PnABcHiBrytpy5bBRx8F1//8Z18Ltskm1ZUxTWcgip2FkJZL2mIqKuEoP+EA+vXXfVA3frzvkBBeWRJgv/18t4gddvAdGm66yd/+4IPlGUDHnYEGP2Fw0KDS7KtjRz+BPT1hcejQ3AF0lOUbkBlAhxfUSgJNIEyMfAuq3jKzgv4HOedexwfJDd0/2zn3PtCCzudVbvTooHH8gAH+m/umm1ZX8CzJlrQMtEo4yk84u3j88X6xjSuvzAyeBw3ynRxGjAgmop10UnD/E08EE6/KSdwZ6Dicc06w/dRTQaeNtPp6f8YxLeoAWhloaUC+AfTuwCgz+8zMxpjZWDMbE+XAwszsbDMbaWYj58yZU6rdJt977wXbSV0xSapbkgLo+vpkZPSkMI1l2TbYwHcM+vBD2H//zPsGDw7qdRct8qUI5SauRVTiNGgQ7Lab3161Cu6+O/P+0aODILJbt+L2/k9Lcg20MtCJkW8AfRAwANgfOBRf23xoVIPK5py73Tk32Dk3uHv37qXabfK9+26wvdNO8Y1DpCFJCqAXLw7q3zt0KL/JstXq3HN9uU2bNr4G95BDfHnGHXf4FezOOCP3v6WZXwwq7cEHSzfmYolrEZW4hScT3n575jyfYcOC7QMOiKYzSZ8+wfb06T6QTwpNIkyMvP7ynHNTgC74oPlQoEvqNomTMtCSdEkKoFW+UZ722QcWLPB/P+PG+UzyDTf49pJN9UUOB9DDhpXfaobVmIEGPzE83QN6yhT4+9+D+6KufwY/mXz99f326tW+VWxSZLexk9jk28buQuBBoEfq8oCZnR/lwKQJ8+bBpEl+u7bWL7kqkjRJCqDVgaN8mTUv07jxxsHZuVWrcre5S7JqzUC3beu/IKWdf76vY4+6fV1YUuuglYFOjHzfkc4EdnLOXemcuxLYGTirsSeY2cPA28CmZjbdzM40s3PN7NzU/eub2XTgYuCK1GP0qZav998PtrfZxp/eFEmaJPWBVgeO6hSeTPjAA/GNozmqcRJh2hVXwHbb+e36ejjhBH9bun3dDjsUv31dWDkE0MpAxyrfNnYGhJsNr07d1iDn3PFN3P8V0CvP/Us2lW9IOUhSBlolHNXp2GPhoov8qfg334QvvkjWynKNqeZJr506+XKN733P94RescK3LkyLqnwjLTyRMCmt7D7+2JczgV/DIbxSp5Rcvhnoe4B3zewqM7sKeAe4K7JRSdM0gVDKQZL6QKuEozr16JF5qv+hh+IbS6GqOQMN/t9uxIjci1odfHC0+05iBjq8JP2hh5Z+aXfJkO8kwuuB0/F9nb8BTnfO3RDlwKQRzikDLeUhqRloBdDV5fjQCdFXX41tGAWr5gx0Wt++PohOTyqE6NrXhSUtgJ45M7OTzM9+Ft9YBGiihMPM1nbOLTSzdYEvU5f0fes65xpcKEUi9MUXMHeu3+7SJeh1KpI0SQ2gVcJRXbbaKthOUkeFplR7Bjpt4EDfRWW//fyXirPPjj77mrQA+uabg4XTdtsNdt453vFIkzXQD+F7Po8CXOh2S13fKKJxSWPC2ecddtBpHEmuJAXQKuGoXhtuGGyXUwBdrW3sctlpJ5g8GT77rDTBY7gX9LRpvoY+rt7xixb5Zc3TlH1OhEYDaOfcIamfZTLjokqofEPKRZICaJVwVK+uXX2nouXL/d/B4sXlEZBWaxu7hnTv7i+l0L69r8GePdu3QJw5E3r3Ls2+s919dzB5cOONff2zxC7fPtAv53OblIgmEEq5SGoArRKO6mKWmYWeOTO+sRRCGeh4JaGMY9Uqv3BQ2sUXaxXVhGg0gDaztqn6525mto6ZrZu61AEbNvZcicjKlfDBB8F1ZaAlyZLUB1olHNWtHMs4lIGOVxIC6CefDPbdtSucemo845A1NFUDfQ5wET5YHkXQ+3khcGuE45KGjB0Ly5b57b59tRKRJFtSM9AKoKtPuBVauWSgNYkwXuEAOo5e0M7Bn/8cXP/xjzPfUyVWTdVA3wjcaGbnO+duLtGYpDGqf5ZykqQ+0CrhqG7lmIFWG7t4hRdTKXUGevVquOyyYNXhNm18AC2JkddKhM65m81sS2Ag0DZ0+/1RDUwaoABaykmSMtAq4ahu4Qx0OQTQK1YEbctat/Yrz0lpxVXCMW+eX7p8xIjgtjPO0NLdCZNXAG1mvwb2xAfQw4CDgDcABdClpgmEUk6SFECrhKO6ldskwuzss1nDj5VoxBFAjx4NRxyRub+DDoJrry3N/iVv+TYQPhrYB/jKOXc6sDXQJrJRSW7ffQeffuq3a2pgu+3iHY9IU5IaQKuEo/qUWwZa9c/xC5dwTJ0K9fXR7u8f/4Bdd80Mnq+4Ap57TiU8CZRvAL3UOVcPrDKztYHZaBGV0ps2zU8qAP8fO9zhQCSJkhJAr16tgKTaldskQtU/x69DB79sOPhymlmzotmPc3DNNfCDH8DSpf62Tp3gqafgt79V27qEyquEAxhpZl2AO/DdOBYD7zX+FCm68H/eDdVFUMpAdgDtXDynorODEa3eWX2ySzjq65P9d6AvfMlQVwdz5/rtL7/M/CJWDMuX+6XJ7w9VxG6yCTzzDGy2WXH3JUWV17uHc+5HzrkFzrmhwH7AqalSDimlcNZEAbSUg1atoG3b4Hq6BWOpqXxD2reHLl389sqVQVCUVFpEJRmirIOeNw/23z8zeN5rL3jnHQXPZSDflQifMbMTzKyDc+5L59yYqAcmOYQz0BtsEN84RAqRhDIOdeAQKK+JhFpEJRmi6gU9YQLsvDO8/npw25lnwgsvwDrrFG8/Epl8z19dD+wOjDOzx8zsaDNr29STpMiUgZZylIRe0OrAIVBeEwmVgU6GKHpBv/Ya7LILTJoU3HbttXDHHWpXWEby7QP9GvCamdUAewNnAXcD1fVJFFf9Zlo4gFYGWspFEjLQKuEQKK+JhMpAJ0OxSzjuvdfXPKd7fLdrBw88AEce2fLXlpLKewaFmbUDjgLOBXYA7otqUIninF+85Iwz/CVOmkQo5ShpAbQy0NWrnFYj1CTCZChWCcfq1XD55XD66UHwvP76Phut4Lks5buQyqPATsALwK3Aq6m2dpVv4sRgwZLWreEPf/B/9HFQCYeUoyQE0KqBFijfEg5loOMTLuGYMqV53VumTYOTT/bBctqgQfD889CnT3HGKSWX71/BPUB/59y5zrl/V03wDL6dzG67+e1Vq+Cuu+IZh3OaRCjlKQkBtEo4BMp3EqEy0PHp1Am6dvXby5fD118X9vzHH4ett84Mng86CN54Q8Fzmcs3gH4d+KWZ3Q5gZgPM7JDohpUw550XbN9+uz8VU2qLFgUTsNq1UxAg5SO84E8SAmhloKtXUxnoVauCxaripgx0cvTrF2yPG5ffcxYv9l01jjkGvvnG39aqFVx5JTz7rN6HKkAhGegVwK6p69OB30UyoiQ66qjgG+jUqTB8eOnHkF2+EedkRpFCJCEDrRIOgcYnEY4e7d9bN9kE5swp7bhyUQY6OQYPDrbfeqvpx3/8MeywA9x9d3Bb374+C3311b4cVMpevgF0f+fcdcBKAOfcUqB6Iri2bTMnEA4dWvoxqHxDylUSAmiVcAhAjx5B/eqcOf6UfNr11/vbJk2C+xIwR15t7JIjXcYJ8OabDT/OObjnHthxRxg/Prj9+OP9F7Tdd49ujFJy+QbQK1JdOByAmfUHljf+lApz9tnB9rBhxV+RqCmaQCjlKmkBtDLQ1at168xJ4F99FWy//Xaw/cEHpRtTQ9TGLjnCAfTbb+cu4/zuOzjtNJ9sW7rU39a+vW9b9+CDwSqYUjHyDaB/je/A0dvMHgReBn4R2aiSaOON/ZKb4L9l3nFHafevHtBSrpKwkIpKOCQtVyu7uXMzF7UYPbq0Y8pFGejkqKsLPncXLvQlGmH19TBkSOaS3AMHwvvvw6mnquSyQuUVQDvnXgKOBE4DHgYGO+dejW5YCXXuucH2nXfCihWl27d6QEu5SloGWiUc1S3XRMJ33818zGefxfe3mqYMdHKYNV7G8e67mV02Tj3Vrx8xcGBpxiexaDSANrPt0hegLzALmAn0Sd1WXQ49NAheZ8+Gp58u3b5VwiHlKmkBtDLQ1S3XRMJ33sl8TH09jBlTujHlokmEyRKuX84OoMOxwEkn+bKNcPchqUhNTQX9SyP3Ofyy3jmZ2d3AIcBs59yWOe434EbgYGAJcJpzLgGFZ41o3RrOOsvPogW47TY49tjS7FuTCKVcJSGAVgmHpOUq4cgOoMGXcey8c2nGlIva2CVLQxlo5+Cpp4Lrxx1XujFJrBrNQDvn9mrk0mDwnHIvcGAj9x8EDEhdzgZuK2TgsTnrLKip8duvvurb1KSX5YySMtBSrpLWB1olHNUtOwO9evWaJRwAH35YujFlq6/PnC+gbGb8tt46SAZMmQLTp/vt8eP9isXg/5322See8UnJNVXC8YvQ9jFZ9/2+sec6514H5jfykMOA+533DtDFzJKfWu3Z05dypJ15pm+yft11sGBBNPt0TpMIpXzFnYFetSrYr5mCkWqXnYEePz4z25sWZwC9ZEmwoEv79kHSRuJTWws77RRcT2ehw+UbBx3k295KVWhqEmH4XMQvs+5rLLucj57AtND16anb1mBmZ5vZSDMbOScJDe6vuSZYWAX8m/All0Dv3nDRRcVvcbdoURAAaBVCKTdxB9DZ9c+aEV/dsicRhss3wqfpx471X77ioPrnZMpVxhEOoA8/vLTjkVg1FUBbA9u5rhcq1/NzrqHqnLvdOTfYOTe4e/fuLdxtEQwc6Gdp//a3sN56we2LF8ONN0L//r42OtdpwebQKoRSzpIUQOvLp2SXcIT7Px96aHD/smX+fT4Oqn9OpuyJhDNm+G4b4OdIHXxwPOOSWDQVQLsGtnNdL9R0oHfoei98h4/y0LUrXHGFzzbfdVdmu5r6enjsMT8BZd994aOPWrYvlW9IOYu7D7Q6cEhY587+TB74v8eXXgru23ln2Hbb4HpcZRzKQCfTzjsHCazRo+GBB4L79twT1lknlmFJPJoKoLc2s4VmtgjYKrWdvj6ohft+FjjFvJ2Bb51zs5p6UuKkl/n++GMYPtwHzGEvv+zfkM8+G77+unn7UA9oKWfhADocGJSKOnBImFlmFnrqVP+zVSsYPDgzgI5rQRUtopJMnTvDoFToU18P114b3KfyjarTVBeOGufc2s65Ts651qnt9PXaxp5rZg8DbwObmtl0MzvTzM41s/RqJMOAz4FJwB3Aj4rw+8THDA480GczRo+GU04JJn6kVy4cMMBPNiy0rk4dOKSc9eoVbE+eXNoFiEAlHLKmXO+jW23lJ5hus01wWxIy0CrhSJZwHfQ33wTb3/9+6cciscp3Ke+COeeOd85t4Jyrdc71cs7d5Zwb6pwbmrrfOed+7Jzr75wb5JwbGdVYSm7rreG++/wklIMOCm5ftMhPNtx/fyhkMqR6QEs5W2cd36kGfPA8blxp968SDsnWM8d89XSHhewSDtfSasVmUAY6ucIBdNrgwb6JgFSVyAJoATbfHIYN86Udm28e3P7KK7D99jAyz+8MykBLudsutHDpByVeL0klHJIt1/toetGUurrgTMU33wQlHqWkDHRyhScSpql8oyopgC6FAw/0EwmvuiqYgDBtmv+PeM89TT9fkwil3MUZQKuEQ7LlykCnA2izzDKOOOqgNYkwufr0WfPvRwF0VVIAXSq1tfDrX8NzzwUf4suX+wmIf2lsxXQ0iVDKX1ICaGWgBdYMgLp0gU02Ca7H3YlDbeySyyyzjGPjjTO7cEnVUABdakOGwPvvwxZbBLf95jcNLweevQqhAmgpR9mdDVavLt2+VcIh2bLfR3fayXfhSIs7gFYGOtnC3baOPlprM1QpBdBxGDDAr36VnnSwcCH85z+5H7twYeYqhAoApBytt16Q9Vu6tLQLVKiEQ7JlZ6DT5RtpcZdwKAOdbKedBhdfDOeeC5dfHvdoJCYKoOPSsWNm25vnn8/9uOzyDX3TlXIVVxmHSjgkW/ZckuwAevPNoU0bvz11KsybV5pxpSkDnWy1tb708rbb9O9TxRRAx+nQQ4Pt557L3S5J5RtSKeIKoFXCIdnatg1aK661Fuy4Y+b9tbWw5ZbB9VJnodXGTiTxFEDHac89feN+gEmTYMKENR+jHtBSKcJ1pXFloFXCIWlDh/qe/HfdBeuuu+b9cZZxqI2dSOIpgI5Tmzb+DTztuefWfIwy0FIpwhnoDz/0S+GWgko4JJf994cXX4STTsp9f5wTCZWBFkk8BdBxyy7jyKYe0FIpevWCbt389sKF8PnnpdmvSjikOcIZ6I8/Lu2+p00Ltrt3L+2+RSQvCqDjdvDBwcTAN9+E+fMz71cPaKkUZvHUQSsDLc2x6abB9sSJpVvSe9EimDHDb9fWBrXaIpIoCqDjtt56wQSW1avhhRcy71cJh1SSUgfQK1bAsmV+u6YG2rePfp9SGbp2hXXW8dtLlmS+F0cp3OJx4419EC0iiaMAOgnCZRzZ7ew0iVAqSakD6Ozss9pASr7MfM/+tIkTS7PfTz8NtjffvDT7FJGCKYBOgkMOIxnzFAAAD2ZJREFUCbaHDw9WJdQqhFJpsgPoqE+Lq3xDWiK8vHeuLklRGD8+2N5ss9LsU0QKpgA6CbbaKliVcMECeOstvx1ehbB9ewUAUv422ihoJTdvXuZkqSiohZ20hDLQItIABdBJYJa7G0d2+YZOP0u5MyttP2h14JCWUAZaRBqgADopwmUcjz0Gf/87PPJIcJvKN6RSlLIOWiUc0hKlDqBXrszMdCuAFkms1nEPQFL22suXaSxZAlOnwrnnZt6vCYRSKeIKoFXCIYUKl3BMnuw7JdXURLe/zz+HVav8dq9eWkRFJMGUgU6Ktm3h2GMbvn+LLUo3FpEoZa9IGCWVcEhLdOoE66/vt1euhClTot1fuP5Z2WeRRFMGOkluugl22gm++MIvqDJ/vp9oNWAAnH9+3KMTKY5NNgnOtsycCccfD5deCltvXfx9qYRDWmrAAPjqK789caKfCBuVcP2zJhCKJJoC6CTp1GnN0g2RSlNTA7vuCv/6l7/+yCP+MmQI/PKXsNtuxduXSjikpTbZBP7zH789YQIccEB0+9IEQpGyoRIOESm9oUNhv/0yb/vnP2H33WGPPfyKnMXoEa0SDmmpUk4kVAs7kbKhAFpESq9/fxgxAkaOhKOOymzR+J//wEEHwfbbwz/+4SduNYdzvgQqTQG0NEepekE7pwy0SBlRAC0i8dl+e3j8cRg3Dk47DVqHqso+/BB+8ANfG/3MM/llpL/9Fu6/H04/Herq4NFHg/tUwiHNUaoM9KxZQclR587B5EURSSQF0CISv802g3vu8a3CLrgA2rUL7vvkEzj8cF83/eqruZ9fX++fv/HGcOqpcO+9vh1k2MCBUY1eKln//sEZkilTYPnyaPaTnX3WwlkiiaYAWkSSo08fuPFGH6hcdllmH9x33vH90nfZBf7wBxgzxmelx4zxddNnnAFz52a+XseOcPDB8NRTmZlEkXy1bev/LsF/Ufv882j2oxZ2ImVFAbSIJE/37nDNNT4jfeGFsNZawX3vvOOD6623ht69fV/pN98M7u/TB373O3j7bd8K8p//9BlskeYqRR20WtiJlBUF0CKSXD16wA03+NrT006DVllvWTNmBJMMa2t9G7xx4+Dyy2Hnnf1tIi1VijpoZaBFykqkAbSZHWhmn5nZJDO7NMf965jZU2Y2xszeM7MtoxyPiJSpvn19jfNXX/lJgscemzkpcK+94KOP4Pe/hw4d4hunVKZSBNDKQIuUlcgWUjGzGuBWYD9gOvC+mT3rnBsXethlwGjn3BFmtlnq8ftENSYRKXPdu8PJJ/vLypXw7ru+LvV739OkK4lO1CUcixb5syngz5pEudqhiBRFlCsR7ghMcs59DmBmjwCHAeEAeiDwBwDn3HgzqzOz9ZxzX0c4LhGpBLW1fuEVkahFnYEOZ58HDMhs5ygiiRRlCUdPYFro+vTUbWEfAUcCmNmOQF+gV4RjEhERKUxdXRDUzpwJixcX9/W1gIpI2YkygM51PjV7JYQ/AuuY2WjgfOBDYNUaL2R2tpmNNLORc+bMKf5IRUREGtK6dWZZxaRJxX19TSAUKTtRBtDTgd6h672AmeEHOOcWOudOd85tA5wCdAe+yH4h59ztzrnBzrnB3bt3j3DIIiIiOYTLOIpdB60JhCJlJ8oA+n1ggJn1M7O1gOOAZ8MPMLMuqfsAfgi87pxbGOGYRERECheeSFjsOmhloEXKTmQzFZxzq8zsJ8CLQA1wt3PuEzM7N3X/UGBz4H4zW42fXHhmVOMRERFptnwmEjoH06bB6NF+Cfpx4/xl0iS/KuYGG8CGG/qfNTUwe7a/hF9PAbRIWYh0qq9zbhgwLOu2oaHtt4EB2c8TERFJlOxWds7BrFkwdiyMGuVbKr73nu9VnsvChX4C4qhRDe+jd+/M5etFJLHUK0dERKQp4Qz0yJHQtSt8803xXr+2Fi5dY70xEUkoBdAiIiJN6dkT2rWDpUv9Ij4NBc+dOsH228OgQTBwoL9suiksW+Yz0LNm+Z/19bDeen65+h49oFevzNU1RSTRFECLiIg0pVUr2GcfeP754La114Ytt4SttoKddoIdd/TBck1N7tfo27c0YxWRyCmAFhERycfDD8Mzz8C66/rAuVcvLSEvUqUUQIuIiOSjY0c48cS4RyEiCRBlH2gRERERkYqjAFpEREREpAAKoEVERERECqAAWkRERESkAAqgRUREREQKoABaRERERKQACqBFRERERApgzrm4x1AQM5sDTIlp992AuTHtu9LoWBaXjmfx6FgWl45n8ehYFpeOZ/FU8rHs65zrnn1j2QXQcTKzkc65wXGPoxLoWBaXjmfx6FgWl45n8ehYFpeOZ/FU47FUCYeIiIiISAEUQIuIiIiIFEABdGFuj3sAFUTHsrh0PItHx7K4dDyLR8eyuHQ8i6fqjqVqoEVERERECqAMtIiIiIhIARRAi4iIiIgUoGoCaDPrbWavmNmnZvaJmV2Yun1dM3vJzCamfq6Tur1r6vGLzeyWrNc63szGmtkYM3vBzLrl2F+tmd2XetynZvbL0vym0SvysfxB6jh+YmbXNbLPX5rZJDP7zMwOiPY3LK1SH08z28/MRqX+NkeZ2d7R/5alEcffZuqxfVKv8bPofrvSi+n/+lZm9nbqcWPNrG20v2VpxPD/vGI/g6BZx7PB9z0z2z51+yQzu8nMrIF9VuTnUKmPZcV8BjnnquICbABsl9ruBEwABgLXAZembr8UuDa13QHYHTgXuCX0Oq2B2UC31PXrgKty7O8E4JHUdnvgS6Au7uOQsGPZFZgKdE9dvw/YJ8f+BgIfAW2AfsBkoCbu41DGx3NbYMPU9pbAjLiPQbkey9DjnwAeA34W9zEo5+OZen8dA2wdel5F/F+P4VhW7GdQM49ng+97wHvALoABw4GDcuyvYj+HYjiWFfEZVDUZaOfcLOfcB6ntRcCnQE/gMPwbEKmfh6ce851z7g1gWdZLWerSIfXNam1gZq5dph7TGmgHrAAWFvWXikkRj+VGwATn3JzU9X8BR+XY5WH4D4LlzrkvgEnAjkX8lWJV6uPpnPvQOZf+m/0EaGtmbYr4K8Umhr9NzOxw4HP8sawoMRzP/YExzrmPUq83zzm3uoi/UmxiOJYV+xkEzTqeOd/3zGwDYG3n3NvOR3T3p5+TpWI/h0p9LCvlM6hqAugwM6vDfwN6F1jPOTcL/B8R0KOx5zrnVgLnAWPxgfNA4K4cD30c+A6Yhc8W/Nk5N784v0FytORY4t+ANjOzutSb/OFA7xyP6wlMC12fnrqt4pToeIYdBXzonFveknEnUSmOpZl1AC4Bri7eyJOpRH+bmwDOzF40sw/M7BfFGn+SlOhYVsVnEDTreIbf93riP1PSGvp8qYrPoRIdy4aeX1aqLoA2s474060XOecK/jZuZrX4AHpbYEP86cZctWU7AqtTj+kH/K+ZbdTccSdRS4+lc+4b/LF8FPgP/hTjqly7yvX0QveXdCU8nun9bQFcC5zTnPEmWQmP5dXAX51zi5s/2uQr4fFsjS9bODH18wgz26eZw06kEh7Liv8MgsKPZ473vXw/Xyr+c6iEx7Kh55eVqgqgU8HvE8CDzrknUzd/nTrtQOrn7CZeZhsA59zk1CmKfwC75njcCcALzrmVzrnZwJtAxawTX6RjiXPuOefcTs65XYDPgIk5HjadzAxLL3KXzZStEh9PzKwX8BRwinNucjF+h6Qo8bHcCbjOzL4ELgIuM7OfFOHXSIwY/q+/5pyb65xbAgwDtivG75EEJT6WFf0ZBIUfzwbe96bjP1PSGvp8qejPoRIfy4r4DKqaADpVr3wX8Klz7vrQXc8Cp6a2TwWeaeKlZgADzax76vp++HqhbFOBvc3rAOwMjG/u+JOkiMcSM+uR+rkO8CPgzhwPexY4LlVj1Q8YgJ+oUBFKfTzNrAvwT+CXzrk3Wzb6ZCn1sXTOfc85V+ecqwNuAH7vnLsl+3HlKob/6y8CW5lZ+1R5wv8A45r/GyRHDMeyYj+DoPDj2dD7Xqo0YZGZ7Zx6zVPI/W9QsZ9DpT6WFfMZ5BIwk7EUF/zpQIcvuRiduhyMn9H8Mv4b/MvAuqHnfAnMBxbjv1kNTN1+Lj5oHgM8B3RN3f594Dep7Y74Wfmf4D8Afh73MUjosXw4dXzGAceFHv/fY5m6fjl+1vNn5JjVW86XUh9P4Ap8beTo0KVH3MehHI9l1r6vovK6cMTxf/0k/Pvmx8B1cR+Dcj2WVPBnUHOOZ2Pve/jM/Mf4z5hbCFZprorPoVIfy8aeX04XLeUtIiIiIlKAqinhEBEREREpBgXQIiIiIiIFUAAtIiIiIlIABdAiIiIiIgVQAC0iIiIiUgAF0CIiMTOzrmY2OnX5ysxmpLYXm9nfItrnRWZ2ShFe5xAzq/ilzEVEwtTGTkQkQczsKmCxc+7PEe6jNfABsJ1zrsHl3rOeU+OcW53jdku91m7Orx4oIlLxlIEWEUkoM9vTzJ5PbV9lZveZ2Qgz+9LMjjSz68xsrJm9kFqKFzPb3sxeM7NRZvZieineLHsDHzjnVplZfzP7ILTPAWY2KrX9pZldaWZvAMeY2QVmNs7MxpjZIwDOZ2FeBQ6J9miIiCSHAmgRkfLRHxgCHAY8ALzinBsELAWGpILom4GjnXPbA3cD1+R4nd2AUQDOucnAt2a2Teq+04F7Q49d5pzb3Tn3CHApsK1zbiv8iqxpI4HvFedXFBFJvtZxD0BERPI23Dm30szGAjXAC6nbxwJ1wKbAlsBLvrKCGmBWjtfZAPg0dP1O4HQzuxj4AbBj6L5HQ9tjgAfN7Gng6dDts4ENm/k7iYiUHQXQIiLlYzmAc67ezFa6YBJLPf793IBPnHO7NPE6S4G2oetPAL8G/g2Mcs7NC933XWh7CLAH8H3gV2a2RaqGum3qNUVEqoJKOEREKsdnQHcz2wXAzGrNbIscj/sU2Dh9xTm3DHgRuA24J9cLm1kroLdz7hXgF0AXoGPq7k2Aj4v1S4iIJJ0CaBGRCuGcWwEcDVxrZh8Bo4Fdczx0OD6THPYg4IARDbx8DfBAqnzkQ+CvzrkFqfv2Av7ZwuGLiJQNtbETEalCZvYU8Avn3MTU9Z8BnZ1zvyrwddYDHnLO7RPBMEVEEkkBtIhIFTKzTYH1nHOvp4Lp/sDezrm5Bb7ODsBK59zoKMYpIpJECqBFRERERAqgGmgRERERkQIogBYRERERKYACaBERERGRAiiAFhEREREpgAJoEREREZEC/D8cmPsnfeSLVAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plot time series using mena values\n",
    "plt.figure(figsize=(12,4))\n",
    "plt.plot(ti, hi,'r',linewidth=3,label='interpolated')\n",
    "plt.xlabel('Time (yrs)')\n",
    "plt.ylabel('Elevation residuals (m)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtkAAAEICAYAAACKx+iJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeZgU1dXH8e9l2PcdFVBQMIKKKIv7K2pUNCoajXEJQU2ixjUa45K4JMZEo4lLEqPRuKLRxGgUFTWo4IIb4IqILAIyLKIssgrMzH3/OF12M/QMPTPVVdXdv8/z1FPdPTVdd+7MdJ8+de69znuPiIiIiIiEp1HcDRARERERKTYKskVEREREQqYgW0REREQkZAqyRURERERCpiBbRERERCRkCrJFREREREIWa5DtnBvunPvEOTfLOXdZLccNcc5VOueOj7J9IiIiIiL10TiuEzvnyoDbgEOAcmCSc26M935aluP+ADyf63N37tzZ9+rVK8TW5mbNmjW0atUq8vMWI/VluNSf4VFfhkv9GR71ZbjUn+Ep5r6cMmXKl977Ltm+FluQDQwFZnnvPwVwzj0CjACmVTvuPOAxYEiuT9yrVy8mT54cVjtzNmHCBIYNGxb5eYuR+jJc6s/wqC/Dpf4Mj/oyXOrP8BRzXzrn5tX4tbhWfEyVfgz33v84dX8ksKf3/tyMY7oD/wQOAu4Gnvbe/6eG5zsDOAOgW7dugx555JE8/wSbW716Na1bt478vMVIfRku9Wd41JfhUn+GR30ZLvVneIq5Lw888MAp3vvB2b4WZybbZXmsesR/C3Cp977SuWyHZ3yj93cCdwIMHjzYx/GJqZg/qUVNfRku9Wd41JfhUn+GR30ZLvVneEq1L+MMssuBnhn3ewALqx0zGHgkFWB3Bo5wzlV475+IpokiIiIiInUXZ5A9CejrnOsNLABOBE7OPMB73zu47Zy7DysXUYAtIiIiIokWW5Dtva9wzp2LzRpSBtzjvf/IOXdW6ut3xNU2EREREZGGiDOTjfd+LDC22mNZg2vv/alRtElEREREpKG04qOIiIiISMgUZItIuObPh1Gj4Lnn4m6JiIhIbBRki0h4li2Dww6DBx6A734XPvgg7haJiIjEQkG2iIRj3To46ij4+GNo1szuf/e7sGJF3C0TERGJnIJsEWm4igo48UR4/XXo2ROmToWBA2H2bPjhD6GqKu4WioiIREpBtog0jPdwzjkwZgx06GC12H36wGOPQfv28NRTcP31cbdSREQkUgqyRaRh/vAHuPNOaN7cAur+/e3x7beHBx+021deCePGxddGERGRiCnIFpH6mzzZAmiAhx+Gfffd9Ovf+Q5cdZWVi4waBZWV0bdRREQkBgqyRaR+1q2zeuuKCrjgAjjmmOzHXXUVdOoEixbBkiXRtlFERCQmCrJFpH5+9SubSWSnneC662o+rqzMBkMCLFgQTdtERERipiBbROpuwgS45RYLoB94AFq0qP347t1tryBbRERKhIJsEamblSvh1FNtVpErroAhQ7b8PQqyRUSkxCjIFpG6ufBCmDcPBg2ykpFcKMgWEZESoyBbRHL36qtwzz22ouMDD0CTJrl9n4JsEREpMQqyRSQ3lZU2iwjAZZel58POhYJsEREpMQqyRSQ3994L774LPXrAJZfU7XuDILu8PPx2iYiIJJCCbBHZsq++gl/+0m7feCO0bFm371cmW0RESoyCbBHZst/+Fr74AvbbD77//bp/f4cOtuz6qlW2iYiIFDkF2SJSu08+gVtvBefS+7pyTtlsEREpKQqyRaR2F11kS6f/6Eewxx71fx4F2SIiUkIUZItIzcaOta1tW7j22oY9V48etleQLSIiJUBBtohkt2GDZbEBrroKunVr2PMpky0iIiVEQbaIZHfbbVaPveOOcN55DX8+BdkiIqWtqgq8j7sVkVGQLSKbW7IEfvMbu33zzdC0acOfU0G2iEjpWr0ahg6FXXaBDz+MuzWRUJAtIpu78kqbG/vww+GII8J5TgXZIiKl69JLYcoUmDYN9t4bnngi7hblnYJsEdnUu+/CXXdB48Zw003hPa+CbBGRktT+nXfgb3+z95Ujj4Q1a+DYY21AfRGXjyjIFikVCxfCccfBm2/WfIz3cMEFtj/vPNhpp/DOv/XWNl/24sU2JaCIiBS/Vav41o032u0rr4QxY+D66+394Mor4fjj4aWXLPAuMgqyRUrFP/8Jjz9u811XVmY/5tFH4dVXoXNnm1EkTE2aQNeuNvBl8eJwn1tERJLpkktosXgx7L47XH65BdeXXmrBdps29r508MHQrp3VbF90ETzzDKxdG3fLG0xBtkipKC+3/bRp8J//bP71tWvhF7+w27/7HbRvH34bVDIiIlI6XnwR7riDqsaN4f77LdkSOPJImDwZzj8fBg2yK6iTJtlg+yOPhI4d4bDD4JZbbKarAiwriTXIds4Nd8594pyb5Zy7LMvXRzjnPnDOveecm+yc2y+OdooUhczA9pprLKOc6Y9/hM8+g4EDLdudDwqyRURKw6pV37yXzB01CnbddfNjdtwRbr3Vgu0VK2DcOLjiChgyxNZq+N//4MILrXRxhx3gnHMKKssdW5DtnCsDbgMOB/oDJznn+lc77EVgN+/9QOB04B/RtlKkiASBbePGm2ez58+3GjmwF7yysvy0QUG2iEhpeOIJmDcPdt+d+SedtOXj27SBb38bfvtbePttKyscPRpOPtmy2nPm2ODJI4+EDh3g//7PAvJx42x6wASKM5M9FJjlvf/Ue78BeAQYkXmA9361999cH2gFFN61ApGkCMpFglUcM7PZl1wC69bBCSfYC1e+KMgWESkNH39s+6OPxtcncdO1K/zgB/DQQ7Z2wxtv2FihIUNg40YbP/S738Ghh1rQPWLElp8zYo1jPHd3YH7G/XJgz+oHOeeOBa4DugLfqenJnHNnAGcAdOvWjQkTJoTZ1pysXr06lvMWI/VluFavXEnVwoU0Al7df3+G3HcfzT/6iI+uuYYNHTuy+yOPUNm0KW8feyzr89jvW61axU7A4nfeYXqB/n71txku9Wd41JfhUn82zM6vvUYXYFpFRXh9eeCBcOCBNF65knZTp9Lugw9o//77tJkxgy9WrmRawn5fcQbZLstjm2Wqvff/Bf7rnPs/4LfAt7M9mff+TuBOgMGDB/thw4aF19IcTZgwgTjOW4zUl+F6/fHHaVRZCZ06sf+RR1oW+6yz2Pk//4FmzQAou+wy9j7xxPw2ZMMGuOEGtqqoYKsC/f3qbxMbgLRypdVQ1jQYyWV7iQdatoQuXb65q/4Mj/oyXOrPBlq+HID+xx7LklWrwu/Lo49O3161iq5ffUXXHj3CPUcDxRlklwM9M+73ABbWdLD3/hXn3A7Ouc7e+y/z3jqRItL0iy/sRlCucdppdpnto4/sfo8eVjKSb8ELoMpFsvvnP632cPjwuFuyud/+Fv71L/jyS1i6tGFznd9/P/zwh+G1TUSSpbISZs602337wjvv5Pd8bdrYljBxBtmTgL7Oud7AAuBE4OTMA5xzfYDZ3nvvnNsDaAosjbylIgWuWfUgu2lT+OUv4ac/tfs33ACtWuW/IZk12d7XnO0sRVOmwCmn2Fyxy5cnq2/Ky+HqqzfNWrdubdM8Zqu1rCm7XVFhiyL94hdwzDHQtm1+2isi8frsM1i/3hYhK+H/89iCbO99hXPuXOB5oAy4x3v/kXPurNTX7wCOA37onNsIrAO+nzEQUkRy1OzL1MWfzEtpp50Gjz1ml+7zXSYSaNvWgvk1a6zcoF27aM5bCG6+2fZffQXLlkGnTvG2J9NDD1ngfPTRcPvt1rZUmVGdeA/77w8TJ8Lvf5+e0UZEisuMGbbfccd42xGzODPZeO/HAmOrPXZHxu0/AH+Iul0ixeabIDvIJIMFSePGRdsQ56wNM2ZYNltBtikvt1KMwPz5yQmyvYcHHrDbP/kJbLNN/Z/LOfswMXSo7c88M5w2ikiyfPKJ7b/1rXjbETOt+ChSArIG2XEJ2hBMKSjw179uWuP82WfxtaW6d96xedW7dLHV1xpqyBAYOdIGwUYxDkBEoqdMNqAgW6QkfFOTnYSR15ore1OrV8Pf/26399jD9kkKskePtv1JJ226JHJD/P730KIF/Oc/tPvgg3CeU0SSQ5lsQEG2SElomsRMtoJsc999NhXePvvA8cfbY0kJsjdutBlPINzZQHr0gEsvBaDPbbelF0USkeIQZLIVZItIsUtkuYiCbJvm6pZb7PZFF8G229rtpATZzz8PX3wB/funs+xhufhi6N6dNjNmpLPlIlL41q6117DGjaFXr7hbEysF2SLFbuVKGq9da5fnO3SIuzUKsjM9/TTMnm1vRMcck7wgOxjwOHJk+FMKtmplZSMAt90W7nOLSHxmzbL9DjuEV2JWoBRkixS7IJjt3j0Zcy8ryE676SbbX3CBzTfdM7U+VxKC7BUrYMwY+5s55ZT8nOOAA2y/aFF+nl9EohfUY5f4oEdQkC1S/DKD7CRQkG1efhleecXmDj/9dHss+CC0cKHVQ8fp0UdtMYmDDkoH/2ELpilcqjXGRIqG6rG/oSBbpNglLcjeaito1AiWLIk/kIzL66/DUUfZ7bPPTq+I1qSJzUPtffwfQoJSkXwuf96qFVVNmsC6dbZlM3Ei9O4N48fnrx0iEh5lsr+hIFuk2AXzUSdh+j6wwTDdulkgWYplAi+/DIceCqtWwQknwDXXbPr1JNRlz5kDr70GLVvCd7+bv/M4x8bgA0ZN2ewxY2DuXLj11vy1Q0TCo0z2NxRkixS7pGWyIR3wx52tjdq4cXD44bas/MiRtlx59YFBSQiy33/f9gccAK1b5/VUWwyyg5lxXnjByldEJLm8VyY7g4JskWKXxCC7FOuyn37aSkTWrYMf/xjuvdey+tUlIcgOAtuttsr7qSpyDbLXrIFXX817e0SkAb780gZNt21rVyxL3BaDbOdcV+fcsc65c5xzpzvnhjrnFJyLFIqklYtAesDbihXxtiMqd9wBI0ZYJvbss22Fx7Ky7McmIcgOAt7g95RHW8xkB6uVAowdm/f2iEgDZGaxkzCbVcxqDJadcwc6554HngEOB7YG+gNXAB86537jnGsbTTNFpN6SmMlu2dL2a9fG2458q6qCSy6Bn/7Ubl9xBfz1rzbwsyZJCLKD7HESguygLaAgWyTpVI+9iSzXKr9xBPAT7/1mr/TOucbAkcAhwGN5apuINNTGjbBkCb5RI1wEl/5zVgpB9rp1MGqUTYXXuDHceSecdtqWvy8JQXYQ8HbunPdTbWzXbtNzVhcE2c2aWZZs9mxb5EJEkkf12JuoMZ3ivf9FtgA79bUK7/0T3nsF2CJJtmgReM+Gjh2z1//GJQiy16yJtx35sngxHHywBdht28Kzz+YWYEM6yJ4/P3/t25IIy0VqrcmuqIDlyy3z/53v2GPPPpv3NolIPSmTvYlcarLbO+fOd87d5Jz7c7BF0TgRaaBUPfb6CDKSdVLMmezJk2HwYHjjDVvEZeJE+Pa3c//+Dh1syfGVK+Grr/LXztoE2eMoMtm1BdnLltm+Y8f0vOIqGRFJLmWyN5HLAMaxQC/gQ2BKxiYiSZeqx1aQHZHRo2G//azf990XJk2CXXap23M4F3/JSFIGPgaDHrt0geHD7fb48cX3dyMSpY8/hr32ghdfDPd5Kypg1iy73bdvuM9doHIJspt77y/y3t/rvb8/2PLeMhFpuFSQvUFBdn5VVMDPf26rI65fD2eeCS+9VP8prIJlzOMKspMy8DEzo77VVjBoEHz9NUyYkPd2iRStu+6Ct96yAdmVleE977x5Ng6oe/e8z69fKHIJskc7537inNvaOdcx2PLeMhFpuKBcpEuXmBtSTTEF2R9+aFmhm26yuvfbb7cp+5o2rf9zxpnJrqqyOmiwMo08q6ht4GP1spWgLlslIyL19/rrtp85E/71r/CeNygVUT32N3IJsjcANwJvkC4VmZzPRolISJJaLtKqle0LOcjeuNGWRB80CKZMscD4pZfgrLMa/txxBtkrVlig3a7d5qtR5kHOmWyAI46w/TPP2MpyIlI369bBO++k7197bXjZ7GDQo+qxv5FLkH0R0Md738t73zu1bZ/vholICJIaZBd6Jvu992DIELj6agu2zzoLpk6F/fcP5/njDLIjHPQIUNGmjd1YvnzzN/ugJjtoy+DBdnvu3HTWTERyN3myvWb16wfbbWf12Y+FNFHc9Om2Vyb7G7kE2R8BBfpOKFLigiBb5SLh2LCBXvfcYwH2++9D7942eOj22yEIFsMQZ5Ad4aBHAF9WBu3bW2a6+gqgQcAf/P2WlaUHQKpkRKTuJk60/QEHwOWX2+1rr7WrVw0xbhzce6/d3mOPhj1XEcklyK4E3nPO/V1T+EmdlZfrsm5cvE9+JruQ5smePBkGDaLX6NE20PG88+CDD+Cgg8I/VwllsoF0QF+9ZCRbW4KSEQXZInUX1GPvuy+ceir06GHjSp58sv7POXEiHHMMbNgA558f3hW9IpBLkP0E8DvgdTSFn9TFk0/aLAnnnBN3S0rT0qU200W7dlS1aBF3azZVaJnsa6+1wY1Tp7K2e3d45RX485/zN4K+Rw/bL1hgAX2UIs5kb3KuXILsffaxfXBpWkRy4306yN5nH1tF9dJL7f4119QvIfbuu/bBd+1aC9pvvtmmIRUghyA7c9o+TeEndfLww7a//XZ46KF421KKUlnsbwK2JCmkIHv2bLjySrucetFFTP7HP/KfqWnWzKasq6y0VTujlPQgOygdWbpUV8lE6mLGDPu/2WorK3UD+PGPYeutbZzJk09avXZFhb3eZfv/8t6+vnatlcwdeqgtnHX88TY1YKNccrelo8becM495Zw7yjm32fBy59z2zrlrnHOn57d5UrAqK61GK3Dmmco8RS01fR/du8fbjmwKKcieN8/2++8Pf/oTVc2bR3PeuEpGklQuUn3gI9jfTvPmNl92Ifz9iCRFZhY7yDY3bw6/+IXdPvZYm3q0SRMb/9CokR1XVmaPNW1qjzVpYjNEDRxorxfDh1sirXHjeH6uBKvtI8dPgP2B6c65Sc65sc65l5xzc4C/A1O89/dE0kopPFOm2JLIvXvDSSdZ7e33vqc3xSgFmeykB9lJz0YGmeStt472vHEF2UnMZFcfuFvT8SJSs2DQ4777bvr4mWfaVKRNmligHATXgaoqy15v3Gj3y8qgRQto29Yy2I891rB1AYpYjR87vPeLgUuAS5xzvYCtgXXADO+9IiWp3fPP2/6ww+CGGyzonjrVBkX84x/xtq1UJDnIbtLEto0bbUvyC3SpBdlJyWSvXWtbs2bpedUDnTvb3/fSpel+EpHaZWayM7VsaQO7q/M+vVVV2dakiUpC6iCn3L73fi4wN+yTO+eGA7cCZcA/vPfXV/v6KUCqKp/VwE+99++H3Q7Jg8wgu00bePRR2HNPuPtu6NrVJqtv1ix9acq59AZ0/PBDmzQ/+DRd7eu0bWtz5paVxfDDFYgk12SDvbB/9ZUFUoUQZG+1VbTnLfVMdnC7c+fNB1IFxwcfCESkdsuW2ZzYzZrlPsVe5nuu3mvrJbYCGudcGXAbcAhQDkxyzo3x3k/LOGwOcID3frlz7nDgTmDP6FsrdfLVV/Dmm3bZKZjebMAA+Mtf4Cc/geuu2+JTDMjlPNtsY6Uop5xitWEa0bypTz+1fRIz2bBpkN2+fdytqVncmez586M9b1Iy2dnqsWs7XkRq9sYbth8yJNlJjSITZ5X6UGCW9/5TAOfcI8AI4Jsg23v/esbxbwIJTcnJJl580QY+7r+/ZZwDP/qRBd4TJ9p8mhs22BRzGzemL0kBeM/SpUvp1KHDJo9tcsynn8KcOfCnP9nWrx/87W8wbFikP2piTZwI48fbi+nuu8PMmXG3aHOFMlf24sW2jzrI7tnT9qWaya6pHhvSgbeCbJHc1FSPLXmVU5DtnGsBbOu9D3Md2+5AZoqmnNqz1D8Cnq3pi865M4AzALp168aECRNCaGLdrF69OpbzJs2O997LNsCnO+7IZ9X7o1cv27Zg9erVtK5tDmLvaTttGt1eeIEu48fT9OOPWXrZZXx4/fU1f0+JcJWVDDrjDFoDc7//febOnJnIv83B3tMamPTyy6wJSlsSaMjs2bQCJpWXs2bChMj6ssny5ewLbPz0UyZG9bvznv/78ksaAS9PnYqfMSPvp1y9ejWTFy1iMLB63jwmp37Wrq+8Qn/g86oqPq728/davZpewJzJk5mXsL/rOCXx/7yQFVN/Dhw7lvbAh23asFTxUXS897VuwFHAJ8Cc1P2BwJgtfV8Oz/s9rA47uD8S+EsNxx4IfAx0yuW5Bw0a5OMwfvz4WM6bKFVV3m+3neWcJ02q99PUqS8//NDO17dvvc+XGFVV3r/xhvdnnOH9qFHeX3ON9w8+aI99+WVuz3HLLdYfvXp5v2aN9z6hf5t77mntfP31uFtSu3btrJ2p/o+sL6uqvG/WzM69cmU051yxws7XunU05/Op/pw3z87bvXv6C7feao+de+7m33Tzzfa1886LrJ2FIJH/5wWsaPpzwwbvW7Sw/5kvvoilCUXTl1kAk30N8WgumexfY6UdE1JB+Xup2UYaqhzomXG/B7Cw+kHOuQHAP4DDvfe6Nph0M2bYvMKdO+c+uKKhdtjB9nPm2DRDhThX5/r18O9/2yqC2UZ5B3r1sqmWBg+22roDDtj05120yBZOAXuuoCQjiQphrux166xuvEkT6Ngx2nM7Z3XZM2daXXb//vk/ZxylIpnny1Yukq0mW+UiIrl77z17Ldtxx2jHWkhOQXaF9/4rF/6gsklAX+dcb2ABcCJwcuYBzrltgceBkd77/F+3lIYLZhU55JDopvlp0cLqV+fPtwA/CLoLxVtvwdFHw5Ildr9jRxsgusMOttpgsE2fDnPn2vbYY3Zs//5Wkz58uN2/+GJYtQqOOsq2JCuEIDuox95qq3gG1vbqZUH2pEnRBNlxDHoE+1to1iy9wEzLlrkNfNTsIiJbpnrs2OQSZE91zp0MlDnn+gLnA69v4Xu2yHtf4Zw7F3gem8LvHu/9R865s1JfvwO4CugE/C0V5Fd47wc39NySR5lT90WpTx8LsmfOLLwge/RoC7D79YOf/xxOPtk+OFRXUWGB9uTJNu/400/DtGlw+OHW39/9Lvzzn7aC1623Rv9z1FUw93GSg+y4ZhYJnHCCrZx6440wcmT+P7jGlcl2zs65cKG1oWVLDXwMQ0WFfXCpq1atNFtTMalpfmzJu1xesc8DdgbWAw8DK4GfhXFy7/1Y7/2O3vsdvPe/Sz12RyrAxnv/Y+99B+/9wNSmADvJ1q+HYGDDoYdGe+6+fW2fxFk0tuSjj2x/0002A0u2ABusLGSXXeDUU206xOnTLfhq184+3Jx5ph33y1/aSptJVwiZ7LiD7B/+0OY5/+gjGDMm/+fLnJs6atVLRmrLqiuTvWUffGCvA23a1H3r0weuvrowX09lc9On2z6qEk75xhaDbO/9Wu/9r7z3Q7z3g1O36/HRWIrea69ZwDRgQPRBSZ8+tp81K9rzhmHqVNvvvHPdvq9ZMysPmTULzj3XFgvYeWf4xS/Cb2M+FFKQHfVCNIGmTdO/z9/9Lv9L0AdBa9SZ7Mxz1iXIViY7u5Urbbnr8nK7stWqVe5bs2Y2Reo111gN7157wX33xf0TSUN8/rnt40oWlLAtBtnOuaecc2OqbaOdcxc455pH0UgpEHGVikDhZrKXLLFgom3b+q/M2LmzZbYXL7ZFgJoXyL9lIcyTHXcmG+DHP7aSicmT4X//y++54ioXyTxnEFzXFmS3bWtXdlavtitokua9/c3MnGkJj2XLrJ9y3dasgRdegFGjoHVrGzNy2mnwvhZbLkiVlfGNtZCcykU+xZY0vyu1rQQ+B3ZM3Rcxr7xi+0MOif7cQZAddybb+7plZoNSkZ13bngNZOfO9qZYKAohkx3XQjSZWra0Wn2wbHY+xflmnJmd9r72rHpQwx0cL2m33QaPPmplH48+WnP5WU3KyuDggy17/fnnNtYD4LnnQm+qRGDpUqiqsv+XJk3ibk3JySXI3t17f7L3/qnU9gNgqPf+HEAFPpIWZJF32SX6c2+/ve2DafziMnq0XXJ9tsZ1kzZV31KRYlAIQXYSMtkAP/2pLT3/6qu25UsSMtlLl9q0iRUVlrFu1mzLxwsAbaZPh4susjv/+IeVezREy5bw/e/b7eBKpRSWoFSkW7d421Gicgmyu6Sm0gO+mVYvSHNsyEurpPAsX26XJVu2jKd+NZjGr6LCpriLS/BG9PjjuR2fmckuNQqyc9e2LZx/vt3OZzY7KQMfc8moa4aRTS1bxs6//jVs3GhjNE44IZznPfhgu3Lw2mvJLu2S7IKpYbt2jbcdJSqXKfx+DrzmnJsNOKA3cLZzrhVwfz4bJwVk9mzb9+kT39RPwTR+s2alB0JGbc4c20+alNvxQZAdR/Y/boUUZMc18DHT+efbnOjPPw+33GJBqXO5bY0aWSZru+3sA0NNUwEmZeBjLkG2ZhjZ1C230Pzzz22Rqj/+Mbzn7dTJnvPtt+Hll+GII8J7bsk/ZbJjtcUg23s/NjU/9k5YkD09Y3aRW/LZOCkgQS10XMEtWF32+PFWthIszhK1IMieOjW9qEZNvC/tcpGkz5NdUWFZIOeS8QbVqZOVjfzxj3DhhfV/nqZN7arP0UfbtJGZklIuUttCNNmOF0swAJx1Vs0lNvV16KEWZP/vfwqyC00QZCuTHYtc158eBPRKHT/AOYf3/oG8tUoKT1KC7My2RG3duvRAucpKG42/9941H79oEaxYAR06JCNTGrWkZ7K/+MI+CHXpkpwBQ1dcYYuLrFhhbct1q6y0hV7mzbOfa/ZsuPlmmwJym23suTMHGxZSuYgy2Wb1atu3aRP+cx96KFx7bf5nt5HwBeUiSUgUlKAtBtnOudHADsB7QGXqYQ8oyJa0ILCNc7XFIMCPaxq/6rXgkybVHmRnloqU4upqSQ+yk1KPnaldO5uusSHWrLErPa+9Bu+8kw6y16616TCfWwwAACAASURBVPCaN6/9Cky+ZAuys632mO14SQfZ+ZhhaK+97Hk//tgy5j17hn8OyQ9lsmOVy8DHwcC+3vuzvffnpbbz890wKTBJymTHFWR/+umm9ydPrv34Ui4VgeTPk52keuwwtWoFe+5pt995J/143HPpauBjw6xaZft8BNlNmsBBB9ntcePCf37JH2WyY5VLkD0VKLJ3GQld5sDHuATT+M2dayPsoxbUYwdL125p8GMpzywCymTHKfgbzQyy46zHBiubcs5KYYKyKw18zF0+y0XASkZAU/kVGg18jFUuQXZnYJpz7vnMVR/z3TApIKtX25tis2b1X7UwDJnT+M2bF/35gyB7xAjL/HzyiS1vXJNSnlkEkh9kJ2EhmnypLciOK5NdVmZzgUP6apQGPuYun5lsSAfZL7xgNf5SGDSFX6xyCbJ/DRwD/B74U8YmYoIs9vbb1zw1WFTiHPwYBNnf+hbstpsNJMsMYjJ5r0x20oPsYs5k9+1rZSPz56dn8ohz+r5AcO7p022vcpHc5TuT3acP9Opl6yHU9LomyeK9arJjtsWIyHv/crYtisZJgUhCPXYgzsGPQZDduzcMHmy3ayoZmT/fMk9dutQ+uKuYKciOT1kZDBxot9991/Zxl4tknnv5ctvnMvBR5SIm35ls5+Cww+y2ZhkpDCtX2mDmVq3SU6ZKpLYYZDvn9nLOTXLOrXbObXDOVTrnarkGLiUnCTOLBOIa/Oh9euDj9tvb4g1Qc5Bd6oMeIfnzZBfrwMdA9ZKRuAc+wuYBfm1tad8+XcNdUZHfdiVdZSWsW4d3Lr8zwwQlIwqyC4MGPcYul2v7fwVOAmYCLYAfpx4TMUnMZEddLrJ8uWUNWrdOr5AGNc8wUur12GCLojRqBBs2JDNIKuZMNqSD7ClTbJ+kTDbY30ZQo51NWRl07Gi3g8x3qUqVilS2aJHf6UAPOsh+L6+/ns6cS3KpVCR2ORXQeu9nAWXe+0rv/b3AsLy2SgpLEmYWCcSVyc4sFXEO+vWzgZhz5mS/nF3q9dhg/ZTUkhHvi3vgI8CgQbYPMtlxD3ysfu6OHS2Qro1KRkwQZOd7fvP27W36x4oKmDAhv+eShlMmO3a5BNlrnXNNgfecczc45y4EVNwjaUnKZO+wgwVvUU/jlxlkAzRuvHmmMJPKRUxSg+wVK6yWsU2b4q1l7NfPFp759FPLBCdp4CPkNlZBM4yYzEx2vh1yiO3Hj8//uaRhlMmOXS5B9sjUcecCa4CewHH5bJQUkHXrbBBf48aw3XZxt8aChh49op/Gr3qQDTUPfqyqspXTQEF2UoPsYq/HBvufHTDAbr/3XjIy2ZlBdi7t0NLqJlW6Udm8ef7PFSxklC15IMmiObJjl8vsIvO8919771d673/jvb8oVT4ikg4ue/WyN+0kiKNkJBj0mBlk1zT4ce5cCyq33jpdU1qqkh5kF2upSCBz8GPSMtm5BNnKZJuoykVg0zKjqqr8n0/qT+UisctldpF9nXPjnHMznHOfBlsUjZMCkKSZRQJxDH4MPmwEq05CzYMfVSqSltQgu9jrsQOZQXbSBj4qyM5dkMmOolykWzfo3t0C+xkz8n8+qT+Vi8Qul3KRu4GbgP2AIRmbSLLqsQNxZLKzlYv06QNt28LChbYFNLNIWlKn8Su1TPbrr8OaNbZSab4WM8mFykXqJ5XJrogikw3pbLZKRpJNmezY5RJkf+W9f9Z7v8R7vzTY8t4yKQxJmlkkEPWqj1VVVgICVjYTaNQoXZedmc3WzCJpSc1kl0qQvcsuFlgHf7+dOuV3Crgt0cDH+okykw3p1zUF2cmmTHbsaiyidc6lUhyMd87dCDwOrA++7r3XuqqSzEx21Ks+Llpkcz136bL5amtDhsBLL8Gjj8L778MLL1jWEBRkQ/KD7GIe+AjQrJkF2sGqj3EOeoT6Z7JLPciOcnYRUCa7UGjgY+xqG6n2p2r3B2fc9sBB4TdHCk4Sg+xgGr85cyx4y/cl1GyDHgNBxufBB9OPNWoEw4en36hKWfC7WbMm3nZUVyqZbLCSkSDIjrMeG2xu+RYtbNaiutRkl3q5SNSZ7OqDHxvltOSGROnrr22BtMaNa1/USfKqxiDbe39glA2RArRhg11mdi57gBmX5s0tuJ00yTLHRx+d3/NlG/QYOPhgm4+4osLml/32t+HAA/WiF0hqJrtUBj6CBdl33223485kgwXO5eUa+FgXUc4uApYZ7dHDfk8zZsBOO0VzXsldUI/dtas+BMWoxp53zl3knPtRlsfPc879LL/NkoIwb55lMbbd1i47J0kQWD/5ZP7PlW3QY6BDB5g2zd6IbrsNjj1WAXampAbZpZbJDsSdyQZbunubbXIL3DTw0USdyQaVjCSdBj0mQm0fb04HRmd5/M7U16TUBYMekzR9X2DECNs/9RRUVub3XLUF2VK7JAbZ69bBV1/ZgMBSmMd8wIB0pisJQfZ999kH+FxmOQl+P8uWlfaczVFnskFBdtJp0GMi1BZke+/9hiwPrgdCGX7unBvunPvEOTfLOXdZlq/v5Jx7wzm33jl3cRjnlBAlsR47sMsuNtPHF1/A22/n91wKsusviUF25qDHOGfaiErLllbSBMkoF3Eu94WtmjSxaTKrquyDUamKeuAjpIPs6usASDIok50ItRbqOOc2++1ke6w+nHNlwG3A4UB/4CTnXP9qhy0Dzgf+GMY5JWRJDrKdS5eMjBmT33PVNvBRapfEebJLqR47sM8+tt9uu3jbUR8qGYm3XOTdd0v7KkJSKZOdCLUF2TcCzzjnDnDOtUltw4CnCCfoHQrM8t5/msqYPwKMyDwgNTf3JGBjCOeTsCU5yIZ0yUg+67LXr4cFC+xy+7bb5u88xSrJmexSCrKvuw7+85/0/0wh0eDH9GI0UQbZweBHrfyYTJq+LxFqDLK99w8AVwLXAHOBOcBvgKu99/eHcO7uwPyM++Wpx6RQJD3I3n9/aNcOPv44f3Nmf/YZeA89e9qla6mbJAfZxT5HdqZOneC443Iv00gSBdnxZLJBJSNJljm7iMSm1ldU7/2zwLN5One2Ykdf7ydz7gzgDIBu3boxYcKE+j5Vva1evTqW88aispL/mz2bRsAr5eVUhfwGF1Zf9hs0iG4vvcSsm2+m/IQTGt6wajpMmsRuwPIOHXg/wb/7pP5tdv70U3YBvpw3j6kJaV+vt9+mFzB33TrmZmlTUvuyUDW0P3eqrGQr4ONXX+XzoPyoxOy1dCnNgZXeR/q3uV2nTvQG5j/5JLN79Aj3yauqcN7jy8pCe0pXUcFWzz5LpzfewGUZEO/8piFI/8pKlqXOX9miBYsOP5xle+5ZEGM1dps+nQ7A+4sXszwBr1el+roZZ9qiHOiZcb8HsLC+T+a9vxOb+YTBgwf7YcOGNahx9TFhwgTiOG/kvIerr7a5n7fdlv87/PDQTxFaX37+Obz0En2mTaNPPn43n3wCQIfdd0/07z6xf5vrbRHZzi1bJqd9//0vAL0GDqRXljYlti8LVIP784knYNw4+nXtSr9S/b1ssDkKmnfuzH5R9sG6dXDPPfT8/HN6Zp73ySfh/vvt/3vjxk23DRs2fyzb40Gd99Ch8L3v2VbfMQPeWznUFVc0qLSlyyuvwO672/Mcc0yy55/eaFW2ux1yiLU5ZqX6uhlnkD0J6Ouc6w0sAE4ETo6xPZIL7+HKK+F3v4OyMvhT9YVBE2b4cLsE/tprdjk57CnKNOixYZJYLrJype3btYu3HZKbUl9a3ft4ZheBTQc/Vlbae8Ljj1tAHMZgSOdsdqi334Zf/MIC7iuugKOOyv053ngDzj8/XdLSpw9cfnnNYy4ystQffPABAwYMsDtTp9r73bvvWmlV//7w97/DfvvV84fLM9VkJ0JsQbb3vsI5dy7wPFAG3OO9/8g5d1bq63c457YCJgNtgarUIjj9vfcr42p3SfMefvUrGyRVVgYPPwzHHx93q2rXrh0MG2YrP44dCyNHhvv8mr6vYZIcZLdtG287JDelvrT6+vV2VbFpU3zU40K6dt105cfPP4eTTrIA+8ILbXXbpk1tvErmlu2x6o+XlVmmfOxYePRRePppC7aPO87Ol0ut8erVlmhZudLGWFx9NfzoRzmPn1nWvLm9f4A9zznnwD33wB/+YIuMXXkljB9f//7Ll8rK9P9Dly7xtqXEbTHIds41A44DemUe772/pqEn996PBcZWe+yOjNuLsTISidvGjZZBuOEGyww/8oi92BWCo4+2IHvMmHCD7FWr4KOP7Ha2JdVly5IYZAfzLSvILgylPvAxlcXOafGefBg82ILeu++Gu+6yso9zzrGsb0Nrl1u2tETO8cfba8Sxx8L//gf33guXXrrl73/pJQuwBw60q5kNrdlv0cJ+tsMOg75901cyk2bpUvug07GjBuTHLJeCoiexqfUqgDUZmxSbjRtttozXX4d//9sy1iefbCvCtWqVDrD/9a/CCbAhPV/2c899UwPcIF98AVddZfWB06ZZXd6OOzb8eUtREufJVrlIYanPPNlr1lg28pxzLFtayIIgu3XreM4flIz86U/2v/O978Gtt4Y/OLBlS7jgArv997/nVo7ybGrehmOOaXiAnWnbbe3nKy//pvY5UVQqkhi5lIv08N4Pz3tLJP+WLbMswNixFkh//bVdVqqosP2KFVYSko1zFkj+8Y91q4dLgu22g912g/fft0/2PXqkt9pqGLO9SaxaZTWHwRvzvvvCr3+djJXyClESM9kqFyksdclkr18Pd94J116bnuJs993hxz/OX/vyLTV9X+xBNsDBB8Po0VbqkQ+HHWav53PmwLhxdr8m3qeD7LAH5zdtCttsY2skLFhgqwsniabvS4xcguzXnXO7eu8/zHtrJD+efBJuvNEGgNT26b9RIxsM0qMHdO9uL2a77mpLlO+8c7iZgKhdfDGce66VAsyY0fDFE4480i5XJnXQS6FIYpCtcpHCkkuQXVlpwd/VV9vVOrDXuAUL4MEHCzvIjrtcZJ99rN65d2+bmadZs/ydq6wMzjjDxgbdcUftQfb06TBvniVABg8Ovy29etnfz9y5yQuylclOjFyC7P2AU51zc4D12PzW3ns/IK8tk3B4D2eeaf90jRvDAQfAEUfAIYdYVrdxY3vhKiuzy+OFuBhFLn7wA9tWrrRLfMGWmvpqM7Vl9Pfbzz54SMM1b277devsA2ASpsRSuUhhyRz46P3mV6DGjbMP2R98YPd33tlmRxo2zILDl1+2wLtQV2yNO5Pdrp31X1lZNP+/p59uH5aeespew2uanzvIYg8fnp929eoFEydakJ00WlI9MXKJqMKfBFmi89ln9g/XsaNdYiv17Fzbtjb1Uv/+cbdEwN78WrSwIHvduvivllRUWFbdufjbIrlp0cKuiKxda1ndIKM7dapN+/bcc3Z/222tTOTkk9PlDEcfbeNPHn44t4F0SRR3JhuiHVy31VY2APLRR22w5dVXZz8uX6UigWDO7iQG2UG5iDLZsdvixzvv/TygPXBUamufekwKwdtv237oUAXYkkxJKhnJrMcugFXdJCXIZo8aZeUL3bpZqdtzz9nv8g9/sIWjRo7ctF74lFNs/+CD0bc5LHFnsuNw1lm2v+su+2Bc3erV8Mor9j986KH5aUNQIjIvgeGQykUSY4tBtnPuAuAhoGtqe9A5d16+GyYheest2++5Z7ztEKlJEoNslYoUlqBk4L//tbEnS5ZYhvvcc2HWLLjkknRpUqbhw+0q39Sp6XKSQpOETHbUDjzQBuIvWADPPLP518ePt1LAoUPzNyg9CLKTnMlWuUjscilU+hGwp/f+Ku/9VcBewE/y2ywJTWYmWySJkhhk66pPYfnrX22Wn7vvhgkTYP58Cz7/8pfaF+No2hROOMFuF2o2uxQz2c6ls9l33LH518emlt/IV6kIJDvIViY7MXIJsh1QmXG/MvWYJF1FBUyZYrcVZEtSJWmubM0sUpj22MNqc08/3QZ39+iR+2C3H/zA9v/8ZzhLgUetFDPZYKVBzZrB889vehUin1P3ZerZ0/bl5dlLVuKkTHZi5PIqdC/wlnPu1865XwNvAnfntVUSjo8+ssBl++01j7MkVxIz2SoXKR377JOeju3ll+NuTd3FvRhNXDp2hFNPtaD6O99JT82Y76n7As2b25S3FRWwcGH+zlNX3iuTnSC5DHy8CTgNWAYsB07z3t+S74ZJCFQqIoUgiUG2Mtmlw7n0AMiHHoq3LfURlIuUWiYb4OabYf/9LZt86KG2Gm+QxT7ssPxPKZjEkpGVK23RpZYtNUNSAtT4F+ica5vadwTmAg8Co4F5qcck6TToUQpBEGSvWRNvO0DlIqUqCLIffdRWwi0kpZrJBhvcOmaMrej7ySdWHvL44/a1fJaKBJIYZAclot27x9sOAWrPZP8ztZ8CTM7YgvuSdMpkSyFIYiZb5SKlpV8/q+teuTL7bBVJVooDHzO1b29TNe6wgwWYEyfa1YnaVoMMSzBXdpKm8bv5ZtsHHxwlVjUG2d77I1P73t777TO23t777aNrotTL6tVWk924Mey+e9ytEalZEoNsZbJLz/e/b/vnn4+3HXVVqgMfM221Ffzvf7aH/E7dlylpmeyPP4ann7Z68bPPjrs1Qm7zZL+Yy2OSMFOm2Ej5AQPskppIUiUpyFa5SOnq18/25eXxtqOuSj2THdh+ewu0DzwQfvnLaM6ZtCD7pptsP2pU7VNXSmRqXFbdOdccaAl0ds51ID1tX1tgmwjaJg2hUhEpFEkKslUuUrq2Sb2tLVgQbzvqKjOTHfz9lqpdd4WXXorufEkKsj//HEaPtlKZCy+MuzWSUmOQDZwJ/AwLqKeQDrJXArfluV3SUBr0KIUiSfNkq1ykdAUDxZI0HVsulMmOz7bb2n7+fKishLKy+Npy2202q8iIEfCtb8XXDtlEbTXZt3rvewMXZ9Ri9/be7+a9/2uEbZT6UCZbCkWSMtkqFyldXbpYkPTllxasFArVZMenRQubi3rjRli0KL52rF1rQTbAz38eXztkM7VlsgHw3v/FObcL0B9onvH4A/lsmDTAokX2ybpNG9hpp7hbI1K7JE3hp3KR0lVWZouLlJfba2hQCpBkVVXp/5vg/0ii1auXlWrMnWsrjcbhvvtg2TJLqu23XzxtkKxyGfh4NfCX1HYgcANwdJ7bJQ0RZLGHDMn/ZPwiDZWkTLbKRUpbUDJSKHXZQYDdurVe6+MSd112ZWV62r6LL7aabEmMXP4rjwcOBhZ7708DdgOa5bVV0jAqFZFCkqQgW+UipS0Y/FgoddmlvBBNUsQ9V/Zf/gKzZkHv3nDssfG0QWq0xXIRYJ33vso5V5FaBXIJoHmyk0yDHqWQJCnIVrlIaSu0TLYGPcYvrkz2+vVwwQXw97/b/V/9ytbFkETJ5Tcy2TnXHrgLm2VkNfB2Xlsl9VdVBZMm2W1lsqUQJCXI3rDBltQuK9Pc8qWq0Kbx06DH+MURZC9cCMcdB2++Cc2awe23w2mnRXd+yVkuAx+DZYPucM49B7T13n+Q32ZJvc2cadm47t3TbxgiSZaUIDuzHlt1jaWp0KbxUyY7flEH2RMnwvHHw+LF0LMnPP44DB4czbmlznIZ+Pikc+5k51wr7/1cBdgJF/yja1YRKRRJmSdbgx6l0MpFlMmOXzBX9mef2ZXkfHrwQTjoIAuwhw2DyZMVYCdcLgMfbwL2A6Y55x51zh2fWg1SkijIwCiLLYUiaZls1WOXrkIb+KhMdvxatbI51jdssOA3H7yHq6+GkSPtPOeeC+PGQdeu+TmfhGaLQbb3/uVUycj2wJ3ACdjgR0miYEJ8BdlSKJIyT7ZmFpHMTLb38bYlF8pkJ0M+S0a+/hpOOQWuucamafzLX2zTIMeCkNPEms65FsBxwFnAEOD+fDZKGiDIwGy9dbztEMlV0jLZCrJLV5s2lplcuzb9oSvJlMlOhnwF2V98AQcfDA8/bL/jp56yLLYUjFxqsv8FfAwcBNwG7OC9Py/fDZN6UrmIFJpgJo+1a+PNHqpcRJwrrMGPymQnQz7myp42zabhff11G+A4cSIccUR4zy+RyCWTfS8WWJ/lvX/Jex9aZb9zbrhz7hPn3Czn3GVZvu6cc39Off0D59weYZ27aAXlIspkS6Fo3BiaNrUAe/36+NqhchGBwprGT4vRJEPYmewXXoB99oE5c2zl5rfeggEDwnluiVQuQfYrwOXOuTsBnHN9nXNHNvTEzrkyLDN+ONAfOMk517/aYYcDfVPbGcDtDT1v3mzYQMu4VnzKpEy2FKIklIyoXESgsDLZKhdJhjCD7LvuguHD7UP/ccfBhAlKmhWwXDPZG4B9UvfLgWtDOPdQYJb3/lPv/QbgEWBEtWNGAA948ybQ3jmXvL+2efNg223Z7eKLoaIivnZ4r0y2FKYkBdkqFylthTSNn8pFkiEIshuSaFu7Fn76UzjjDKishEsvhX//O/3aKAUpl+GpO3jvv++cOwnAe7/OuVBWaugOzM+4Xw5UXwc82zHdgUXVn8w5dwaW7aZbt25MmDAhhCbmyHuGNm1Ky/nzmXr99Xy5337RnTtDk6++Yt+NG9nYujUTg6XVC9Tq1auj/R0WuaT359BGjWgJvDV+POt69oylDX2nTaM7MPPzz1lQS18lvS8LTdL6s/uaNfQFFkyaxMwEtSubnefMoQswde5cvpwwIXF9Wehy7c+ydevYH6icM4dXx4+v82JWrWfNot+119Jq3jyqmjRhxgUXsHj4cHjllfo1PIFK9W8zlyB7Q2p2EQ/gnNsBCKNwMttfYfVRT7kcYw96fyc2xSCDBw/2w4YNa1Dj6uzCC+Gii9jl1VfhiiuiPXfgA1snqMm22xL5zx+yCRMmFPzPkCSJ789OnaC8nD133RUGDoynDXffDUDfQYPoW0tfJb4vC0zi+vPLL+G22+gOdE9Su7JpbktW7LL33jBsWPL6ssDVqT+7dKHsiy8Y1rt3OrO9JVVVcMstcPnlNv91v340+uc/2WngQIptOblS/dvMpVzkauA5oKdz7iHgReCSEM5dDmSmrHoA1YvgcjkmGUaNorJpU/jf/2D27HjaoDmypVAlYa5slYsI1D7wcfp0+POf87+yX65Uk50cQ4bY/vXXczv+iy9stpCf/9wC7J/+1FZwjCvJIHmRy2I044DvAqcCDwODvfcTQjj3JKCvc663c64pcCIwptoxY4AfpmYZ2Qv4ynu/WalIInTsyBcHHmi377wznjZojmwpVEmoydbsIgK1D3w891y44AJ4+ulo21QT1WQnR1AmOnHilo995RULpp9/3q7iPfkk/O1vqr8uQjUG2c65PYIN2A6rg14IbBvGVHre+wrgXOB5bB7uf3vvP3LOneWcOyt12FjgU2AWcBdwdkPPm08Ljj7abtxzTzxTkWlmESlUSQiyNbuIQDpJsXixDUALVFTAm2/a7UmTom9XNspkJ8e++9q+tiC7qgquuw4OPNDer/fbD957D4LYQYpObTXZf6rlax5bnKZBvPdjsUA687E7Mm574JyGnicqq/r1s0+n770Hjz0GJ58cbQM0s4gUqiQF2SoXKW1Nm0KXLnY5f8mS9OvpRx+ly5nefTe+9mVSJjs5hgyBJk3gww/tqlj115GNG+GYY2BsKuS57DL47W+1PHqRqzGT7b0/sJatwQF2UXLO6qoAbo9hSm9lsqVQJSnIViZbsk3jF2SxwRIpSaDFaJKjRQvYYw/LVmf+rQTGjLEAu2NH2193nQLsElBbucglGbe/V+1rv89nowraySdbVuG112Dq1GjPrYGPUqhatbK9arIlCbINfswMnBYssEx3nDZssK1JE2jWLN62iKmtLvuJJ2x/2WVw+OHRtUliVdvAxxMzbl9e7WvD89CW4tC6NYwcabfvuKP2Y8OmgY9SqOLOZK9fnw5YUtOiSQnLNvgxCLLbt7d93CUjymInT0112Rs3pgfLHnNMtG2SWNUWZLsabme7L5nOSo3bfOABWLEimnNqtUcpZHEH2ZmlIqGstSUFrXome/lym76veXM44QR7LO6SEQ16TJ59Ugtjv/WWBdaBV16xWKB/f+jbN562SSxqC7J9Dbez3ZdMu+5qn2hXrYLevW2i+WzTQYVp6VL7p+7QwWrDRApJ3PNkq1REMlXPZL/9tu0HDYKhQ+12UjLZGvSYHN26QZ8+9jr2/vvpx4NSEWWxS05tQfZuzrmVzrlVwIDU7eD+rhG1r3DddZcF2itWwPXX2wpQo0Z9sypj6FQqIoUsKZlszSwisHkmOygV2Wsv2H13ux13kK1MdjJVr8v2XkF2CattdpEy731b730b733j1O3gfpMoG1mQ+vWzwY9vvAHHH2/zrT7wAOy2GxxyCDz3nP3zhUWDHqWQJSXIViZbYPPZRTKD7J13tlkhZsyId4VSZbKTqXpd9jvvQHm5/U0NGhRfuyQWuSyrLg2x117w6KMwcyacf77NovDCCza6eNddLfDOXPCgvjR9nxSyuINslYtIpsxykaoqq7EFez1v1sxqa73P35XJXCiTnUyZQXZmFnvECGikkKvU6Dcele23h1tvhfnzrXxkm21scYNRo2wS+5dfbtjzq1xEClkQZAfZuaipXEQydepkM80sX26B9PLl9prdo4d9PQklI8pkJ9O3vmVzYS9cCHPnqlSkxCnIjlqHDnDppTBnDtx7r71ov/suDBtmZSWfflq/51W5iBSyXr1sP2NGPOdXuYhkatQo/Vr6+OO232uv9NeDIDvOGUaUyU6mRo3S2ewHHrD1Mtq1gwMOiLddEgsF2XFp2hROPRU++QSuucYyeY89ZrXcf/tb3eu1lcmWQtavn02PNnt2dNNeZlK5iFQXBNmPPWb7zCB74EDbK5Mt2QRB9o032v4737H3fCk5CrLj1rIlXHmlZfBGjrQFMc45garY0QAAEkdJREFUB04/Hb7+OvfnUSZbClnjxjBggN2OIzuochGpLqjLnjbN9tmC7A8/3HQ+5ChpMZrkCoLsYGCsSkVKloLspOje3S4tPfSQzXN9332w//7w2We5fb8GPkqh22MP27/zTvTnVrmIVBcE2QBlZZvODNGunY2zWb/erkbGQeUiyTV4cDpz3bQpDNci2aVKQXbSnHyyTfvXuzdMnmwv7G+8Ufv3aLVHKQZxBtkqF5HqMhMWu+2WHpwbiLtkROUiydW8uQXaAN/+tn5HJUxBdhLttpsF2IceCl9+CWefXfvxmas9Nm8eTRtFwpaETLbKRSSQmcnOLBUJxD3DiDLZyTZihO1HjYq3HRIrBdlJ1bEjjBlj82q/917tZSMqFZFisMsuVps9fXr0i3yoXESqy3w9rS3IjmuGEWWyk+3nP7f1MU44Ie6WSIwUZCdZs2aWzQZ45pmaj9PMIlIMmjWz1fS8h/ffj/bcKheR6raUyc4sFwlz9d5cKZOdbGVl0KdP3K2QmCnITrojj7T9U0/VfIxmFpFiEVfJiMpFpLqePW1Rmh12yB4sbbMNdOliU07Omxd9+5TJFkk8BdlJ953v2P6ll2q+hK5MthSLIMiOus5V5SJSXYsWMGUKvPIKOLf5152Lt2REmWyRxFOQnXTdusHQoTZV1AsvZD9GmWwpFnFksr1XuYhkt912tb+uxjnDiDLZIomnILsQHHWU7Z9+OvvXNfBRisVuu1mGcOpU+2AZha+/hooKqwlv1iyac0pxyFyUJkobNsCSJbaEd8eO0Z5bRHKmILsQZAbZVVWbf13lIlIsWrWCnXayoHfq1GjOqVIRqa8dd7T9zJnRnnfWLKistPUUNG2rSGIpyC4EAwbYIJzFi61GsDqVi0gxibpkRKUiUl99+9p+1qzsCZB8mT7d9v36RXdOEakzBdmFwLn0LCPVS0a02qMUm6iDbM0sIvXVtq2Nm/n6aygvj+68H39s+512iu6cIlJnCrILRU1T+Wm1Ryk2cQXZymRLfQQlIzNmRHfOIJOtIFsk0RRkF4qDDoKWLW0U+4IF6cc16FGKTTCY7P337QNkvqlcRBoijrrsIJOtchGRRFOQXSiaN4dDDrHbmSUjGvQoxaZ9e1sAZP36dMYun1QuIg0R1GVHlcn2XplskQKhILuQZKvL1qBHKUZRloyoXEQaIupykfJyW5isa1dN3yeScLEE2c65js65cc65mal9hxqOu8c5t8Q5F9FcXgkXrP749NP24tqnD1x5pT2mIFuKSbCSXhSLfKhcRBoiyGRHVS6iLLZIwYgrk30Z8KL3vi/wYup+NvcBw6NqVOJtvTX84Ad2e/lymD07XZ+9667xtUskbHFkslUuIvWxww42A9Snn0YzhkD12CIFI64gewRwf+r2/cAx2Q7y3r8CLIuqUQVh9GhbqOPLL+3y5BtvwJtvwoknxt0ykfAEmewpU+D5560ONV9ULiIN0aIFbLutLQ4zd27+z6dMtkjBaBzTebt57xcBeO8XOee6NvQJnXNnAGcAdOvWjQkTJjT0Kets9erVsZwXgFdeiee8eRJrXxahQuzPgbvuSvsPP4Thw1nVty+fnXIKX+y3H5SVhXqefjNn0g2YVl7Okhz6qBD7MsmKoT8HdO5Mx3nz+OA//2HZ3nvn9Vy7vfkmHYAPNmxgWbV+K4a+TBL1Z3hKti+993nZgBeAqVm2EcCKascur+V5egFT63LuQYMG+TiMHz8+lvMWI/VluAqyP7/6yvvrrvO+a1fvLZft/Y47en/33d6vXx/eeY44wp77qadyOrwg+zLBiqI/zz7b/oZuuin/59pqKzvX3Lmbfako+jJB1J/hKea+BCb7GuLRvJWLeO+/7b3fJcv2JPC5c25rgNR+Sb7aISIFqm1buOwyuwR/222w3XZWIvWjH1kd7K232iwLDaVyEWmoqAY/rlgBixfbmgk9e+b3XCLSYHHVZI8BRqVujwKejKkdIpJ0LVrA2WdbAPPAA9C/v01j9rOfWeB9/fWwdm3dnnPtWnjhBfjlL+HDD+0xBdlSX1FN4xfUY3/rW9BIM/CKJF1c/6XXA4c452YCh6Tu45zbxjk3NjjIOfcw8AbwLedcuXPuR7G0VkTi16QJjBxpQfETT8DQobB0KVx+uU1nefvtW57d4b33YPhw6NDBFne67jqbwq9rV9h++2h+Dik+Ua36qEGPIgUlliDbe7/Ue3+w975var8s9fhC7/0RGced5L3f2nvfxHvfw3t/dxztFZEEadQIRoywWXWefx4GDbJFmc4+26Y1u/deWFZtUqKVKy3zPWiQfc/GjTZN4MUXw9ixNh2mMtlSX716QePG8NlnsG5d/s6j6ftECkpcs4uIiDSMc3DooZaRfvxx+NWv4JNP4PTTLRDfe29bwKlLF7jqKgvEGzWyYPtXv4LOneP+CaRYNG5sV0JmzIBZs/K3boEy2SIFRUG2iBQ25+C44yy7PXo0PPQQvPwyTJxoW2DvveFvf4OBA+NrqxSvHXe0IHvmzPwF2UEmW0G2SEFQkC0ixaFxYzjtNNtWrrSBjc88Y9nt4HENFpN8CWYYydfgx/XrbVXJRo3S5xKRRFOQLSLFp21b+O53bROJQr4HP86aZatK7rADNG+en3OISKiU1hEREWmofGeyg3psDXoUKRgKskVERBoq35lsDXoUKTgKskVERBqqe3dbOOnzz23u9bBp+j6RgqMgW0REpKEaNbJFkSA/2WxlskUKjoJsERGRMOSrZKSqSkG2SAHS7CIiIiJhqOvgx+XLYdo0277+Grbe2rZttoH27W3l0iVLLGhfswa6doWOHfPXfhEJlYJsERGRMNSWyV64EN5+27ZJk2DqVFi8uG7Pv/PODW+jiERGQbaIiEgYgiD73XfhwQctkJ46Fd57DxYs2Pz4li1tIOPOO0OrVrBokW0LF8KKFdCpk2Wvu3aFbt3g7LOj/XlEpEEUZIuIiIQhKBeZNg1Gjtz0a23bwtChtg0ZAgMHwrbbahVSkSKmIFtERCQMXbrAD34AkyfDLrvYtvPOsOuuFoAroBYpKQqyRUREwuAcjB4ddytEJCH0sVpEREREJGQKskVEREREQqYgW0REREQkZAqyRURERERCpiBbRERERCRkCrJFREREREKmIFtEREREJGQKskVEREREQua893G3IXTOuS+AeTGcujPwZQznLUbqy3CpP8OjvgyX+jM86stwqT/DU8x9uZ33vku2LxRlkB0X59xk7/3guNtRDNSX4VJ/hkd9GS71Z3jUl+FSf4anVPtS5SIiIiIiIiFTkC0iIiIiEjIF2eG6M+4GFBH1ZbjUn+FRX4ZL/Rke9WW41J/hKcm+VE22iIiIiEjIlMkWEREREQmZgmwRERERkZApyE5xzvV0zo13zn3snPvIOXdB6vGOzrlxzrmZqX2H1OOdUsevds79tdpzneSc+9A594Fz7jnnXOcs52vinLs/ddzHzrnLo/lJoxFyf34/1ZcfOeduqOWclzvnZjnnPnHOHZbfnzA6Ufelc+4Q59yU1N/mFOfcQfn/KaMTx99m6thtU89xcf5+umjF9H8+wDn3Ruq4D51zzfP7U0Ynhv/1on0fqkdf1vi655wblHp8lnPuz845V8M5i/I9CKLvz6J5H/Lea7O69K2BPVK32wAzgP7ADcBlqccvA/6Qut0K2A84C/hrxvM0BpYAnVP3bwB+neV8JwOPpG63BOYCveLuhwT2ZyfgM6BL6v79wMFZztcfeB9oBvQGZgNlcfdDgfbl7sA2qdu7AAvi7oNC7s+M4x8DHgUujrsPCrUvU6+vHwC7ZXxfUfyfx9SfRfs+VI++rPF1D3gb2BtwwLPA4VnOV7TvQTH1Z1G8DymTneK9X+S9fyd1exXwMdAdGIG9QJHaH5M6Zo33/jXg62pP5VJbq9Sns7bAwmynTB3TGGgBbABWhvpDxSjE/twemOG9/yJ1/wXguCynHIG9Waz33s8BZgFDQ/yRYhN1X3rv3/XeB3+zHwHNnXPNQvyRYhXD3ybOuWOAT7H+LBox9OWhwAfe+/dTz7fUe18Z4o8Uqxj6s2jfh+rRl1lf95xzWwNtvfdveIv4Hgi+p5qifQ+C6PuzWN6HFGRn4ZzrhX2Kegvo5r1fBPZHBnSt7Xu99xuBnwIfYsF1f+DuLIf+B1gDLMIyDn/03i8L5ydIlob0J/ZCtZNzrlfqjeAYoGeW47oD8zPul6ceKyoR9WWm44B3vffrG9LupIqiP51zrYBLgd+E1/Lkiehvc0fAO+eed86945y7JKz2J01E/VkS70P16MvM173u2PtJoKb3lpJ4D4LI+rOm7y8oCrKrcc61xi7r/sx7X+dP9M65JliQvTuwDXZpM1ud21CgMnVMb+Dnzrnt69vupGpof3rvl2P9+S/gVexyZkW2U2X79rqeL8ki7MvgfDsDfwDOrE97ky7C/vwNcLP3fnX9W5tsEfZlY6w84pTU/ljn3MH1bHZiRdifRf8+VNe+zPK6l+t7S9G/B0Gk/VnT9xcUBdkZUgHyY8BD3vvHUw9/nrq8QWq/ZAtPMxDAez87dSnk38A+WY47GXjOe7/Re78EmAgMDuHHSIyQ+hPv/VPe+z2993sDnwAzsxxWzqaZmh5kL9MpSBH3Jc65HsB/gR9672eH8TMkScT9uSdwg3NuLvAz4JfOuXND+DESIYb/85e9919679cCY4E9wvg5kiLi/izq96G69mUNr3vl2PtJoKb3lqJ+D4LI+7Mo3ocUZKek6qfvBj723t+U8aUxwKjU7VHAk1t4qgVAf+dcl9T9Q7Dapeo+Aw5yphWwFzC9vu1PmhD7E+dc19S+A3A28I8sh40BTkzVfPUG+mKDKwpe1H3pnGsPPANc7r2f2LDWJ0/U/em9399738t73wu4Bfi99/6v1Y8rRDH8nz8PDHDOtUyVQRwATKv/T5AsMfRn0b4P1bUva3rdS5VArHLO7ZV6zh+Svf+L9j0Iou/Ponkf8gkYfZmEDbv06LHyjvdS2xHYKO0XsSzAi0DHjO+ZCywDVmOfzvqnHj8LC6w/AJ4COqUePxq4JnW7NTbTwEfYm8Qv4u6DBPfnw6k+mgacmHH8N/2Zuv8rbET3J2QZrVyoW9R9CVyB1Wm+l7F1jbsfCrU/q5371xTX7CJx/J//AHvdnArcEHcfFHJ/UsTvQ3Xty9pe97Ds/lTs/eWvpFfLLon3oDj68//bu58Xq8o4juPvT2PgQiiIkIJAmNKFFZUYqBWly5GEMNwF/gMRIdLGcuOiCFoEuYl+gIYtIheKP4K0aNmMk1NJiOAmDCkoMLSm/La4Z+iQp0I7l+bOvF9w4XnOc873nHMXlw+H557nn44fpY/LqkuSJEk9c7qIJEmS1DNDtiRJktQzQ7YkSZLUM0O2JEmS1DNDtiRJktQzQ7YkzXNJbksy3Xy+S/Jt076U5I0hnfO5JM/0UGdzkgW9pLwkdfEVfpI0QpLsBi5V1atDPMcSYAp4qKq6luPuOmasqn7v2J6m1oYarNIoSYuCT7IlaUQleTzJoaa9O8m7SY4nOZ/kqSSvJJlJcrRZEpkka5J8kmQyybG5JZH/YiMwVVW/JRlPMtU65z1JJpv2+SQvJvkMeDrJs0m+TnI6yQGAGjzJOQlsHu63IUnziyFbkhaOcWAC2ALsA05U1X3AZWCiCdqvA1urag3wFrCno84GYBKgqs4BPyV5oBnbDrzT2vdKVT1SVQeAF4AHq+p+BivfzvkceLSfW5Sk0bDk/74ASVJvjlTVbJIZYAw42myfAVYAq4B7gY8GszgYAy501LkDONPqvwlsT/I8sA14uDX2fqt9Gtif5CBwsLX9InDnDd6TJI0kQ7YkLRy/AFTV1SSz9eefbq4y+L0P8FVVrfuXOpeBpa3+B8BLwMfAZFX90Br7udWeAB4DngR2JVndzOle2tSUpEXD6SKStHh8A9yeZB1AkpuTrO7Y7wxw91ynqq4Ax4C9wNtdhZPcBNxVVSeAncCtwLJmeCXwZV83IUmjwJAtSYtEVf0KbAVeTvIFMA2s79j1CIMn0m37gQKO/035MWBfM1XlFPBaVf3YjD0BHP6Ply9JI8VX+EmSrpHkQ2BnVZ1t+juAW6pq13XWWQ68V1WbhnCZkjRvGbIlSddIsgpYXlWfNoF7HNhYVd9fZ521wGxVTQ/jOiVpvjJkS5IkST1zTrYkSZLUM0O2JEmS1DNDtiRJktQzQ7YkSZLUM0O2JEmS1LM/AF5fILL34XVyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#conver time series ot elevation  change with respect to the start date\n",
    "plt.figure(figsize=(12,4))\n",
    "plt.plot(ti, hi - hi[0],'r',linewidth=2)\n",
    "plt.ylabel('Elevation Change (m)')\n",
    "plt.xlabel('Time (yrs)')\n",
    "plt.grid()\n",
    "plt.savefig('time_series_3km_v2.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/scipy/ndimage/filters.py:1447: RuntimeWarning: Mean of empty slice\n",
      "  cval, origins, extra_arguments, extra_keywords)\n"
     ]
    }
   ],
   "source": [
    "#estimate dh/dt (assuming it's linear)\n",
    "nt,nx,ny = Zi.shape\n",
    "rate = np.zeros((nx,ny))*np.nan\n",
    "for i in range(nx):\n",
    "    for j in range(ny):\n",
    "        dh = Zi[:,i,j]\n",
    "        if np.any(np.isnan(dh)):continue\n",
    "        rate[i,j] = np.polyfit(ti,Zi[:,i,j],1)[0]\n",
    "Zf = generic_filter(rate.copy(), np.nanmean, 2) #smooth the data - filtering parameter can be messed with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trend: 0.0 (m/yr)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAIBCAYAAABHgbZKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de7SddX3n8ff35E4uBAi3BBCIwQoKFlloq61YLxMoll7sFJ3ebEd7w+q0a7W2dI0dV1enHaeuVnFK01YrrdV2VqWgRhBdVmpbKwmDSEAhUpQQbgFJIPdzznf+2Dt6jOdk7+ec3z5n//K8X2vtlbP3fp7vfs4m+s3n+f2e3xOZiSRJao+RuT4ASZI0u2z+kiS1jM1fkqSWsflLktQyNn9JklrG5i9JUsvY/CVJmoGIOD0iPhMR90TEloh4yyTbRES8OyK2RsSdEXHhXBzrIfPn8sMlSToKjAK/npm3R8RyYHNE3JKZd0/Y5lJgXffxIuBPu3/OCZO/JEkzkJkPZ+bt3Z+fBu4B1hy22RXAddnxeWBlRJw6y4f6TTZ/SZIKiYgzge8G/v2wt9YAD054vo3v/AfCrPG0vySpaseddHEePLBzoJ+xe+e9W4B9E17akJkbJm4TEcuAfwDempm7DisRk5Sds/X1bf6SpKodPLCTF3zfht4bzsC/fOySfZl50VTvR8QCOo3/g5n5kUk22QacPuH5acD2skfZP5u/JKlqQRAjczeKHREB/CVwT2a+a4rNbgSuiogP05notzMzH56tYzyczV+SpJl5CfBTwJci4o7ua78NnAGQmdcCG4HLgK3AHuANc3Cc32TzlyTVLSBGJhtSnx2Z+TkmH9OfuE0CvzI7R9Sbs/0lSWoZk78kqXoRZtkm/LYkSWoZk78kqXIxp2P+NTL5S5LUMiZ/SVLdgjm9zr9GfluSJLWMyV+SVLUARhzzb8TmL0mqXHipX0N+W5IktYzJX5JUtzle3rdGJn9JklrG5C9Jqp6X+jXjtyVJUsuY/CVJVQvH/Bsz+UuS1DImf0lS5YIRr/NvxG9LkqSWMflLkqrnmH8zJn9JklrG5C9Jqpuz/Rsz+UuS1DImf0lS1cK7+jXmtyVJUsuY/CVJ1XPMvxmTvyRJLWPylyTVLWDEu/o14rclSVLLmPwlSdWLcMy/CZu/JKlqQTjhryFP+0uS1DImf0lS3VzetzGTvyRJLWPylyRVb8TlfRvx25IkqWVM/pKkqgWO+Tdl8pckqWVM/pKkuoXX+Tdl8pckqWVM/pKk6rm8bzMmf0mSWsbkL0mq3ohj/o2Y/CVJahmTvySpahEQI2bZJvy2JElqGZO/JKly4Wz/hkz+kiS1jMlfklQ9V/hrxuQvSVLLmPwlSXULr/NvyuYvSapa4PK+TXnaX5Kkljli8l+/fn3u2LFjto5FknSU2bx5882ZuX7QnzPXE/4i4n3A5cBjmfm8Sd6/BLgB+I/uSx/JzHfM3hF+uyM2/x07drBp06bZOhZJ0lEmIlbN9THMkr8CrgGuO8I2/5yZl8/O4RyZY/6SpMoFI3M85p+Zt0bEmXN6EA045i9J0uz4noj4YkR8IiLOm8sDMflLkuoWszLmvyoiJo6Db8jMDQ32vx14VmY+ExGXAf8IrCt6hA3Y/CVJ6m1HZl403Z0zc9eEnzdGxP+JiFWZOSez6m3+kqSqBXM/27+XiDgFeDQzMyIupjPs/sRcHY/NX5KkGYqIDwGX0Bke2Aa8HVgAkJnXAq8FfikiRoG9wJWZmXN0uDZ/SVL95nqFv8x8XY/3r6FzKeBQcLa/JEktY/KXJNXNG/s0ZvKXJKllTP6SpMrFnI/518bkL0lSy5j8JUlVq+E6/2Fj8pckqWVM/pKkugUY/JuZleZ//1e/WrTeM+PLi9YDeGJf2Zp3fX1x0Xp3fvHJovWefc5xResBvOr8bxStdzLbi9ZbMLq/aD2AYLxsvfGxovU2jb6waL29B+YVrQfw0mWbi9cs6alFJxevOUbZ73Eh5f9un712bfGaGh4mf0lS9Rzzb8Yxf0mSWsbkL0mqnsm/GZu/JKlqEcGIi/w04ml/SZJaxuQvSaqep/2bMflLktQyJn9JUvUc8m/G5C9JUsuY/CVJVYuAEcf8GzH5S5LUMiZ/SVL1wkH/Rkz+kiS1jMlfklQ9V/hrxuQvSVLLmPwlSdULo2wjfl2SJLWMyV+SVLWIcLZ/QyZ/SZJaxuQvSaqeK/w1Y/KXJKllZiX5Z+GxmIUjB4rWA1gwb7xovfGy5Xj2OccVrXfGqeX/lbxndHHRektGnypab9HuJ4rWA3hi1XOK1ltyYFfReqWdsHRf8ZoL95X9ne9delHRel/bsaJoPYBnrSz7O6898EDRet2qA6g5OA75N2PylySpZRzzlyTVLSAc82/E5i9JqloA9v5mPO0vSVLLmPwlSdVzkZ9mTP6SJLWMyV+SVLeAEaNsI35dkiS1jMlfklS1wDH/pkz+kiS1jMlfklS9MMo24tclSVLLmPwlSZULRhzzb8TkL0lSy5j8JUl1C2/p25TJX5KkljH5S5Kq1rnOf66Poi4mf0mSWsbkL0mq3siI0b8Jk78kSS0zK8k/MovWWxAHitYDOGb+vqL1zl69qGi9JQvGitZbOG+0aD2A/WPzitZb/I2HitbLe+8qWg9gyfeeWrTegQXHFK13DAeL1kvKp6tHjn1O0XoPPrGiaL1BBMqzDt5TtN4xD99XtB4Az39p+ZoD5Jh/MyZ/SZJaxuYvSapaROcMzSAfvY8h3hcRj0XEpKcYo+PdEbE1Iu6MiAtLfw9N2PwlSdWLkRjoow9/Baw/wvuXAuu6jzcBfzrjX3oGbP6SJM1QZt4KPHmETa4ArsuOzwMrI6LspKEGvNRPklS9WZjwtyoiNk14viEzNzTYfw3w4ITn27qvPVzi4Jqy+UuS1NuOzLxoBvtP9s+TspfCNWDzlyRVr4I1frYBp094fhqwfY6OxTF/SZJmwY3AT3dn/b8Y2JmZc3LKH0z+kqTKxRDc0jciPgRcQmduwDbg7cACgMy8FtgIXAZsBfYAb5ibI+2w+UuSNEOZ+boe7yfwK7N0OD3Z/CVJ1RtxELsRvy5JklrG5C9Jqt5cj/nXxuQvSVLLmPwlSdULo38jJn9JklrG5C9Jqlr0edtdfYvJX5KkljH5S5Kq55B/MyZ/SZJaxuQvSaqeyb+ZKpv/grH9xWseO39n0XpLVu4rWm8xe4vWmz9+oGg9gEfGVxett3/lqUXrLV5yf9F6AAv3P1203v4FS4vWO3Ph14vWO2b/U0XrAcTesrc0v+i4Z4rWW/XY3UXrAYx/9t+L1hstWq3r1YMoqmFRZfOXJOmQwLX9m7L5S5LqNgS39K2N/1aSJKllTP6SpOq5yE8zJn9JklrG5C9Jqp5j/s2Y/CVJahmTvySpaoHJvymTvyRJLWPylyRVz9n+zZj8JUlqGZO/JKlurvDXmMlfkqSWMflLkqrWme1f9u6QRzuTvyRJLWPylyRVz9n+zZj8JUlqGZO/JKl6zvZvxuQvSVLLmPwlSdUz+TczK80/C/9XmZejResBHLPnqaL19i1cXrTeM/NWFq23ZOzpovUAli3YXbTeQ0ufU7TeSRcuLVoPYNeiE4rWG2de0XprHru9aL0d7/+rovUAjl17WtF6Jz33uUXrPf2FTUXrASz7rnVF6+0/78VF6+noZ/KXJFUtAka8zr8Rx/wlSWoZk78kqXqO+Tdj85ckVc/m34yn/SVJahmTvySpaoHL+zZl8pckqWVM/pKk6gVe6teEyV+SpJYx+UuSquds/2ZM/pIktYzJX5JUt3C2f1Mmf0mSWsbkL0mqWpCEN/ZpxOQvSVLLmPwlSdVztn8zJn9JklrG5C9Jqt6IK/w1YvKXJKllbP6SpOpFDPbR+/NjfUR8JSK2RsTbJnn/kojYGRF3dB//fRDfQ79m5bR/ZNnTMQdHFhWtB7B8//ai9XJkXtF6Xxs9s2i9PfOXFK0HMHqw7F+nPaMLi9bbwbFF6wEcfKbsv5/XLttWtN68R79etN43HniiaD2A8dHxovVO+r4fKFpv9CffUrQewJZ4VtF6259eUbQewA8Vr3j0ioh5wHuBVwHbgNsi4sbMvPuwTf85My+f9QOchGP+kqSqBcz1df4XA1sz836AiPgwcAVwePMfGp72lySpt1URsWnC400T3lsDPDjh+bbua4f7noj4YkR8IiLOG+jR9mDylyTVbXbW9t+RmRdNfQTf4fBTEbcDz8rMZyLiMuAfgXUlD7AJk78kqXpBDvTRwzbg9AnPTwO+bSJZZu7KzGe6P28EFkTEqpLfQRMmf0mSZuY2YF1EnAU8BFwJvH7iBhFxCvBoZmZEXEwnfE9rBm1EnNb9jO8DVgN7gbuAjwOfyMyes2ht/pKk6s3l8r6ZORoRVwE3A/OA92Xmloj4xe771wKvBX4pIkbpNOsrM5tfChcR76czn+BjwB8CjwGLgXOA9cDVEfG2zLz1SHVs/pIkzVD3VP7Gw167dsLP1wDXFPioP8rMuyZ5/S7gIxGxEDijVxHH/CVJVTt0qd8gH8MiM++KiHkR8TdTvH8gM7f2qmPzlySpIpk5BpzYTfnT4ml/SVLlso039nkA+JeIuBHYfejFzHxXPzvb/CVJqs/27mMEWN50Z5u/JKl6cznbf478wxQT//rimL8kSfW5NiK+EBG/HBErm+5s85ckVa8ts/0PycyXAj9JZ2XBTRHxtxHx6n73t/lLklShzLwX+B3gN4GXAX8SEV+OiB/tta9j/pKkqgX0s/7+USUizgfeAPwgcAvwmsy8PSJWA/8GfORI+9v8JUmqzzXAXwC/nZl7D72Ymdsj4nd67WzzlyTVbXZu6TsUImID8AngBzPz6cm2ycy/7lXHMX9JkurxPuACYGNEfDoifjMiLmhaxOQvSapeW8b8M/PzwOeB342IE4BXA7/enQNwO3BTZv59rzqz0vzHC59gGMbLLg634MCeovVWLpv07M60jUTP2z03NjYyr2i95+y5rWi97cc+t2g9gGuuX1q03n++9NSi9VYvmPbS35M657+8qmg9gJETTixa76YFVxSt9+BXy58gPf/Mvb03auDEpbt7b9TYigHUVEmZ+QTwoe6DiHghndv69mTylyRVr4ZQWFJ3YZ+fBs5kQi/PzF/tZ3+bvyRJ9dlI5/T/l4DGp3Jt/pKkqgXZmjH/CRZn5q9Nd2ebvySpem077Q/8dUS8EfgYsP/Qi5n5ZD872/wlSarPAeCdwNXwzdMeCZzdz842f0lS9Vq4aM2vAc/OzB3T2bmF35ckSdXbAkz7mnKTvySpei0c8x8D7oiIz/DtY/5e6idJ0lHqH7uPabH5S5Kq1sZb+mbmB2ayv2P+kiRVIiI+GhGviYgFk7x3dkS8IyJ+rlcdk78kqXLZpjH/N9KZ6f/HEfEk8DiwGDgL2Apck5k39Cpi85ckqRKZ+QjwG8BvRMSZwKnAXuDezOx79r/NX5JUvbaN+QNk5gPAA9PZ1zF/SZJaxuQvSapbtPI6/xkx+UuSVKGIWBIRz5nOvjZ/SVL1Dt3Wd1CPYRMRrwHuAG7qPn9BRNzY7/42f0mS6vO7wMXAUwCZeQdwZr87O+YvSapaACNDmM4HbDQzd0bEtHaeleafhU8wjOSBovUADixaUbTeSI4VrbdnbEnRerd/bWXRegDPPnVf0XqrN/7fovXOev7zi9YDOOmUXyha70tfK/s/yd2n/nzRemd/16NF6wGMZtnf+aZ/2FW03tLlC4vWA3jZObuL1jtxbHvReh0vGEBNFXRXRLwemBcR64BfBf6135097S9Jqlxnhb9BPobQm4Hz6NzR70PALuCt/e7saX9JkirTXc3v6u6jMZu/JKl6wzgjf5Ai4qPwHb/0TmAT8GeZecRxWE/7S5Kq17ZL/YD7gWeAP+8+dgGPAud0nx+RyV+SpPp8d2Z+/4TnH42IWzPz+yNiS6+dbf6SpKoF7TvtD5wYEWdk5tcBIuIMYFX3vZ6XxNn8JUmqz68Dn4uIr9L5989ZwC9HxFLgA712tvlLkqo3pJfjDUxmbuxe3/9ddJr/lydM8vvjXvvb/CVJqtML6SzpOx84PyLIzOv62dHmL0mq3NDOyB+YiPhrYC2dm/scWlI2AZu/JElHqYuAczNzWv/qsflLkqrXtuQP3AWcAjw8nZ1t/pIk1WcVcHdEfIHO+v4AZOYP9bOzzV+SVL0WJv/fncnONn9JkiqTmZ+dyf6u7S9JqtqhFf7atLZ/RLw4Im6LiGci4kBEjEXErn73t/lLklSfa4DXAfcBS4D/2n2tLzZ/SVLlkmB8oI9eImJ9RHwlIrZGxNsmeT8i4t3d9++MiAtn/FtnbgXmZeZYZr4fuKTffR3zlyRpBiJiHvBe4FXANuC2iLgxM++esNmlwLru40XAn3b/nK49EbEQuCMi/hedS/6W9ruzyV+SVL2IHOijh4uBrZl5f2YeAD4MXHHYNlcA12XH54GVEXHqDH7ln6LTw68CdgOnAz/W784mf0mSelsVEZsmPN+QmRu6P68BHpzw3ja+M9VPts0aprlIT2Z+rfvjPuB/NN2/yuZ/IBYXrxmLVhatt3BsX++NGhgdL3uS5ulneo9hNXVwrOwxLnrFpUXrfeO4M4vWA3jFeM/bZjey92DZ/y7/fEcUrffFJTMJKpN78XPLfocxsr/3Rg2MjZWf6b1wpOzv/M2V3VssprfKbRM7MvOiqT5+ktcOP6B+tulbRLyEzrX+z2JCL8/Ms/vZv8rmL0nSENlG57T7IacB26exTRN/Cfw3YDPT+OefzV+SVL05vhb/NmBdRJwFPARcCbz+sG1uBK6KiA/TGRLYmZnTOuXftTMzPzHdnW3+kiTNQGaORsRVwM3APOB9mbklIn6x+/61wEbgMmArsAd4w3Q+a8Ilgp+JiHcCH+Hb1/a/vZ86Nn9JUtWCJLL8PKYmMnMjnQY/8bVrJ/ycwK8U+Kg/Ouz5xHkICfxAP0Vs/pKk6g3jEryDkJkvL1HH6/wlSapERPxaRPz8JK+/OSLe2m8dk78kqXpzfdp/Fv0cMNnSwBvoTDz8436KmPwlSapHdlcRPPzF/Uy+lsCkTP6SpOq1ZcwfICJOzsxHD3+tSQ2TvyRJ9Xgn8PGIeFlELO8+LgE+CvzvfouY/CVJdcu5v9RvtmTmdRHxOPAO4Hl0Lu/bAry9yaI/Nn9JkirSbfLTXt0PbP6SpMoF7RrzL8Exf0mSWsbkL0mq3izc0veoYvOXJKkyEbEI+DHgTCb08sx8Rz/72/wlSZVrz2z/CW4AdgKbmXBXv37Z/CVJqs9pmbl+ujs74U+SVL0gB/oYQv8aEc+f7s4mf0mS6vNS4Gcj4j/onPYPOuv+n9/PzrPS/NetfVbRenfc93jRegAxUna8aMXuR3tv1MC6+VuK1lt9/vFF6wGs2P1I0Xr3HP+yovUe272saD2AtceW/Z1X7yn733nkwkuK1vvkv40VrQewc9/CovVecMGSovXWrBotWg/ghH0PFa03kuX/u9SmhWP+l85kZ0/7S5JUmcz8GrASeE33sbL7Wl9s/pKkqgWd6/wH+Rg2EfEW4IPASd3H30TEm/vd3zF/SZLq8/PAizJzN0BE/CHwb8B7+tnZ5i9JqlwStG7MP4CJkz3Guq/1xeYvSapbAkN4an7A3g/8e0Rc333+w8Bf9ruzzV+SpMpk5rsi4p/oXPIXwBsy8//1u7/NX5JUvbZc6hcRKzJzV0QcDzzQfRx67/jMfLKfOjZ/SZLq8bfA5XTW9J841hHd52f3U8TmL0mq3NAuwVtcZl7e/fOsmdTxOn9JkioTEZ/u57WpmPwlSdVr0Zj/YuAYYFVEHMe3Lu9bAazut47NX5KkevwC8FY6jX4z32r+u4D39lvE5i9Jql9LrvPPzD8B/iQi3pyZfa3mNxmbvyRJlcnM90TE84BzgcUTXr+un/1t/pKkqgXZmjH/QyLi7cAldJr/Rjq3+P0c0Ffzd7a/JEn1eS3wCuCRzHwDcAGwqN+dTf6SpLolQ3nb3QHbm5njETEaESuAx+hzgR+w+UuSVKNNEbES+HM6s/6fAb7Q7842f0lS/Vo25p+Zv9z98dqIuAlYkZl39ru/Y/6SJFUmIm6IiNdHxNLMfKBJ44dKk/8L1p1YvOb9X91VtF6Mjxatt3fRSUXrPTx6StF6AAeWLe69UQPL2F203vHH9nWzq0bGCv9PaNeKNUXrnTd+T9F6uy88v2g9gNHx6L1RA+edVvbvzepFjxStB7Dg6b1F6+1ZcnzRevVp32x/4F3ATwD/MyK+APwd8LHM3NfPzlU2f0mS2iwzPwt8NiLmAT8AvBF4H51lfnuy+UuSqteWu/pNFBFLgNfQOQNwIfCBfve1+UuSVJmI+DvgRcBNdNb0/6fM/sc+bP6SpPq1b8z//cDrM3NsOjs721+SVLWgs8jPIB9D6FbgtyJiA0BErIuIy/vd2eYvSVJ93g8cAL63+3wb8Hv97uxpf0lS3TLbeNp/bWb+RES8DiAz90ZE39fNmvwlSarPge5s/wSIiLXA/n53NvlLkuo3nOPyg/R2OjP9T4+IDwIvAX62351t/pIkVSYzb4mI24EX05nz+JbM3NHv/jZ/SVL12rK8b0RceNhLD3f/PCMizsjM2/upY/OXJKkef3SE95LOUr892fwlSZXL1oz5Z+bLS9Rxtr8kSZWIiN+Y8POPH/be7/dbx+YvSape5PhAHzM6tojjI+KWiLiv++dxU2z3QER8KSLuiIhNU5S7csLPv3XYe+v7PSabvyRJg/U24NOZuQ74dPf5VF6emS/IzIumeD+m+Hmy51Oy+UuS6pZ0Vvgb5GNmruBbt9v9APDDM/xtJ/t5sudTcsKfJEm9rTrsVPyGzNzQ574nZ+bDAJn5cEScNMV2CXwyIhL4synqXxARu+ik/CXdn+k+X9zn8dj8JUm1m5U77+04wql4IuJTwCmTvHV1g894SWZu7/7j4JaI+HJm3jpxg8yc16DelGz+kiTNUGa+cqr3IuLRiDi1m/pPBR6bosb27p+PRcT1wMV0bt1bnM2/6+y1a4vWe/C+vu+v0Jd9I0uL1hsZK1oOgB0HJp3AOm2nLnikaL2FY/uK1gN4Yt7JRevNn3ewaL1VT32laL3zVh1btB7AvvElResdzLL/t7a//zOpfVuw96mi9UaXlv17WKXxoV7h70bgZ4A/6P55w+EbRMRSYCQzn+7+/GrgHYM6ICf8SZI0WH8AvCoi7gNe1X1ORKyOiI3dbU4GPhcRXwS+AHw8M28a1AGZ/CVJ9RviFf4y8wngFZO8vh24rPvz/cAFs3VMJn9JklrG5C9JqltmiWvxW8XmL0mq3ixc6ndU8bS/JEktY/KXJNXP0/6NmPwlSWoZk78kqXJO+GvK5C9JUsuY/CVJVQuc7d+UyV+SpJYx+UuS6pYM+419ho7JX5KkljH5S5Iql0N9Y59hZPKXJKllTP6SpPp5nX8jJn9JklrG5C9Jqp9j/o2Y/CVJahmTvySpbple59+QzX9ATl93btF6j96zuWi93fOPKVoP4PE9K4rWeyROLlpv2fy9ResBrN731aL1lux8uGi9PceuLlrvhF1fL1oPYNGTDxatN7r8hKL1Hj7+eUXrATy56pyi9ZYcfLpoPR39bP6SpPo5278Rx/wlSWoZk78kqX7O9m/E5i9JqpsT/hrztL8kSS1j8pck1c/T/o2Y/CVJahmTvySpfl7q14jJX5KkljH5S5LqN+6YfxMmf0mSWsbkL0mqW6Zj/g2Z/CVJahmTvySpfq7w14jJX5KkljH5S5Lq5wp/jZj8JUlqGZO/JKlyzvZvyuQvSVLLmPwlSXVLXOGvIZt/JU5+7guL1vvSnfuL1gM4fdnjRestG3uqaL0n86Si9QB2LSlbc8+ilUXrHRxZVLTeyY9+sWg9gCg8UWvBNx4tWu+k+YuL1gPYuXxN0XrLHr2vaD0Azr24fE0NDZu/JKlqCaRj/o045i9JUsuY/CVJlUvH/Bsy+UuS1DImf0lS3RKv82/I5i9Jql56Y59GPO0vSVLLmPwlSZVLb+zTkMlfkqSWMflLkuqWgGP+jZj8JUlqGZu/JKl+mYN9zEBE/HhEbImI8Yi46AjbrY+Ir0TE1oh424w+tAebvyRJg3UX8KPArVNtEBHzgPcClwLnAq+LiHMHdUCO+UuSKpdDfZ1/Zt4DEBFH2uxiYGtm3t/d9sPAFcDdgzgmk78kSb2tiohNEx5vKlx/DfDghOfbuq8NhMlfklS3ZDZu7LMjM480Xv8p4JRJ3ro6M2/oo/5kpwUG9kvZ/CVJmqHMfOUMS2wDTp/w/DRg+wxrTsnmL0mqXtZ/Y5/bgHURcRbwEHAl8PpBfZhj/pIkDVBE/EhEbAO+B/h4RNzcfX11RGwEyMxR4CrgZuAe4O8zc8ugjsnkL0mq3+DH/KctM68Hrp/k9e3AZROebwQ2zsYx2fxb6qxlDxWvuea+z5QtuPMbRcstet73F60H8OD8tUXrLRw5WLTeMewuWm//8pOK1gNYPOQ3ZFn4zBPFax4/dqBswfAkrpqx+UuS6pYJ9Y/5zyr/uShJUsuY/CVJ1cshHvMfRiZ/SZJaxuQvSarfEK/tP4xs/pKkqmUmOeRXjQwbT/tLktQyJn9JUv087d+IyV+SpJYx+UuSquelfs2Y/CVJahmTvySpbi7v25jJX5KkljH5S5Kq55h/MyZ/SZJaxuQvSapeep1/IyZ/SZJaxuQvSapbAo75N2LylySpZUz+kqSqJUl6nX8jNv+WWnv22eWLFq65+89/p2i9pY/fX7QeQKwu+zvvG1tUtN5x8XjRemPzFhatBxAHD5QtWLgJjC4/oWg9gH1Ljita7/GV5xetB3Bu8YoaJjZ/SVL1vM6/Gcf8JUlqGZO/JKluCXidfyMmf0mSWsbkL0mqXDrm35DNX5JUt3R536Y87S9JUsuY/CVJ1cv0tH8TJn9JklrG5C9Jqlx6qV9DJn9JklrG5C9Jqlqmy/s2ZfKXJKllTP6SpOqZ/Jsx+UuS1DImf0lS9VzhrxmTvyRJLWPylyTVLb2xT1Mmf0mSWsbkL0mqnmP+zZj8JUlqGZO/hsPvsR8AAASNSURBVNbSN/7eXB9CT9891wfQ03FzfQC9nXvxXB/BrFteuN6JhevVxhX+mjP5S5LUMiZ/SVL90uTfhMlfkqSWMflLkiqXzvZvyOYvSaqbE/4a87S/JEktY/KXJFXO0/5NmfwlSWoZm78kqXo5ngN9zERE/HhEbImI8Yi46AjbPRARX4qIOyJi04w+tAdP+0uSNFh3AT8K/Fkf2748M3cM+Hhs/pKkug378r6ZeQ9ARMz1oXyTp/0lSRoOCXwyIjZHxJsG+UEmf0lS9WZhtv+qw8bhN2TmhkNPIuJTwCmT7Hd1Zt7Q52e8JDO3R8RJwC0R8eXMvHUGxzwlm78kSb3tyMwpJ+tl5itn+gGZub3752MRcT1wMWDzlyTpO+TMZ+TPtYhYCoxk5tPdn18NvGNQn3fE5r958+bNwzRBQZJUnYHPXB92EfEjwHuAE4GPR8QdmfmfImI18BeZeRlwMnB9t+fOB/42M28a1DEdsfkf6RSHJEnDYnxseJN/Zl4PXD/J69uBy7o/3w9cMFvH5Gx/SZJaxjF/SVLVOtf5u7Z/EyZ/SZJaxuQvSape7bP9Z5vJX5KkljH5S5IqV/91/rPN5C9JUsuY/CVJdRvyu/oNI5O/JEktY/KXJFXP6/ybsflLkqqWnvZvzNP+kiS1jMlfklS5HOob+wwjk78kSS1j8pck1c0x/8ZM/pIktYzJX5JUPS/1a8bkL0lSy5j8JUlV8zr/5kz+kiS1jMlfklQ5r/NvyuQvSVLLmPwlSXVzzL8xk78kSS1j8pckVc/r/Jsx+UuS1DImf0lS1TIhne3fiMlfkqSWMflLkirndf5NmfwlSWoZk78kqW5e59+YzV+SVLUET/s35Gl/SZJaxuQvSapbQo65yE8TJn9JklrG5C9Jqlw64a8hk78kSS1j8pck1S2d7d+UyV+SpJYx+UuSqpZ4Y5+mTP6SJLWMyV+SVLeE8VGTfxMmf0mSWsbkL0mqW0IeNPk3YfKXJKllTP6SpKplpmP+DZn8JUlqGZO/JKlujvk3ZvKXJKllTP6SpLp5nX9jJn9JklrG5C9JqltCHhyf66OoislfklS1Q5f6DfIxExHxzoj4ckTcGRHXR8TKKbZbHxFfiYitEfG2GX1oDzZ/SZIG6xbgeZl5PnAv8FuHbxAR84D3ApcC5wKvi4hzB3VAnvaXJNVtyC/1y8xPTnj6eeC1k2x2MbA1M+8HiIgPA1cAdw/imEz+kiTNnp8DPjHJ62uAByc839Z9bSBM/pKkqm1l/82Xj967asAfszgiNk14viEzNxx6EhGfAk6ZZL+rM/OG7jZXA6PAByfZLiZ5bWCnM2z+kqSqZeb6ITiGVx7p/Yj4GeBy4BWZOVlT3wacPuH5acD2ckf47TztL0nSAEXEeuA3gR/KzD1TbHYbsC4izoqIhcCVwI2DOiabvyRJg3UNsBy4JSLuiIhrASJidURsBMjMUeAq4GbgHuDvM3PLoA4oJj/7IEmSjlYmf0mSWsbmL0lSy9j8JUlqGZu/JEktY/OXJKllbP6SJLWMzV+SpJax+UuS1DL/H0FW2NBe9CQ2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x648 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plot spatial distribution of dh/dt\n",
    "fig = plt.figure(figsize=(9,9))\n",
    "ax = plt.axes(projection=ccrs.UTM(7)) #or UTM(7)\n",
    "im = plt.pcolormesh(Xi, Yi, Zf, transform=ccrs.UTM(7), cmap='coolwarm_r',vmin=-1.5,vmax=0.5)\n",
    "\n",
    "plt.colorbar(label='Elevation Change (m/yr)')\n",
    "plt.clim([-2, 2])\n",
    "print('Trend:',np.around(np.nanmean(Zf),1),'(m/yr)')\n",
    "plt.savefig('gridded_ts_3km_filteted.png')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function generic_filter in module scipy.ndimage.filters:\n",
      "\n",
      "generic_filter(input, function, size=None, footprint=None, output=None, mode='reflect', cval=0.0, origin=0, extra_arguments=(), extra_keywords=None)\n",
      "    Calculate a multi-dimensional filter using the given function.\n",
      "    \n",
      "    At each element the provided function is called. The input values\n",
      "    within the filter footprint at that element are passed to the function\n",
      "    as a 1D array of double values.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    input : array_like\n",
      "        The input array.\n",
      "    function : {callable, scipy.LowLevelCallable}\n",
      "        Function to apply at each element.\n",
      "    size : scalar or tuple, optional\n",
      "        See footprint, below. Ignored if footprint is given.\n",
      "    footprint : array, optional\n",
      "        Either `size` or `footprint` must be defined.  `size` gives\n",
      "        the shape that is taken from the input array, at every element\n",
      "        position, to define the input to the filter function.\n",
      "        `footprint` is a boolean array that specifies (implicitly) a\n",
      "        shape, but also which of the elements within this shape will get\n",
      "        passed to the filter function.  Thus ``size=(n,m)`` is equivalent\n",
      "        to ``footprint=np.ones((n,m))``.  We adjust `size` to the number\n",
      "        of dimensions of the input array, so that, if the input array is\n",
      "        shape (10,10,10), and `size` is 2, then the actual size used is\n",
      "        (2,2,2). When `footprint` is given, `size` is ignored.\n",
      "    output : array or dtype, optional\n",
      "        The array in which to place the output, or the dtype of the\n",
      "        returned array. By default an array of the same dtype as input\n",
      "        will be created.\n",
      "    mode : str or sequence, optional\n",
      "        The `mode` parameter determines how the input array is extended\n",
      "        when the filter overlaps a border. By passing a sequence of modes\n",
      "        with length equal to the number of dimensions of the input array,\n",
      "        different modes can be specified along each axis. Default value is\n",
      "        'reflect'. The valid values and their behavior is as follows:\n",
      "    \n",
      "        'reflect' (`d c b a | a b c d | d c b a`)\n",
      "            The input is extended by reflecting about the edge of the last\n",
      "            pixel.\n",
      "    \n",
      "        'constant' (`k k k k | a b c d | k k k k`)\n",
      "            The input is extended by filling all values beyond the edge with\n",
      "            the same constant value, defined by the `cval` parameter.\n",
      "    \n",
      "        'nearest' (`a a a a | a b c d | d d d d`)\n",
      "            The input is extended by replicating the last pixel.\n",
      "    \n",
      "        'mirror' (`d c b | a b c d | c b a`)\n",
      "            The input is extended by reflecting about the center of the last\n",
      "            pixel.\n",
      "    \n",
      "        'wrap' (`a b c d | a b c d | a b c d`)\n",
      "            The input is extended by wrapping around to the opposite edge.\n",
      "    cval : scalar, optional\n",
      "        Value to fill past edges of input if `mode` is 'constant'. Default\n",
      "        is 0.0.\n",
      "    origin : int or sequence, optional\n",
      "        Controls the placement of the filter on the input array's pixels.\n",
      "        A value of 0 (the default) centers the filter over the pixel, with\n",
      "        positive values shifting the filter to the left, and negative ones\n",
      "        to the right. By passing a sequence of origins with length equal to\n",
      "        the number of dimensions of the input array, different shifts can\n",
      "        be specified along each axis.\n",
      "    extra_arguments : sequence, optional\n",
      "        Sequence of extra positional arguments to pass to passed function.\n",
      "    extra_keywords : dict, optional\n",
      "        dict of extra keyword arguments to pass to passed function.\n",
      "    \n",
      "    Notes\n",
      "    -----\n",
      "    This function also accepts low-level callback functions with one of\n",
      "    the following signatures and wrapped in `scipy.LowLevelCallable`:\n",
      "    \n",
      "    .. code:: c\n",
      "    \n",
      "       int callback(double *buffer, npy_intp filter_size,\n",
      "                    double *return_value, void *user_data)\n",
      "       int callback(double *buffer, intptr_t filter_size,\n",
      "                    double *return_value, void *user_data)\n",
      "    \n",
      "    The calling function iterates over the elements of the input and\n",
      "    output arrays, calling the callback function at each element. The\n",
      "    elements within the footprint of the filter at the current element are\n",
      "    passed through the ``buffer`` parameter, and the number of elements\n",
      "    within the footprint through ``filter_size``. The calculated value is\n",
      "    returned in ``return_value``. ``user_data`` is the data pointer provided\n",
      "    to `scipy.LowLevelCallable` as-is.\n",
      "    \n",
      "    The callback function must return an integer error status that is zero\n",
      "    if something went wrong and one otherwise. If an error occurs, you should\n",
      "    normally set the python error status with an informative message\n",
      "    before returning, otherwise a default error message is set by the\n",
      "    calling function.\n",
      "    \n",
      "    In addition, some other low-level function pointer specifications\n",
      "    are accepted, but these are for backward compatibility only and should\n",
      "    not be used in new code.\n",
      "\n"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
